Piecewise linear neuron modeling

ABSTRACT

Methods and apparatus for piecewise linear neuron modeling and implementing one or more artificial neurons in an artificial nervous system based on one or more linearized neuron models. One example method (for implementing a combination of a plurality of neuron models in a system of neural processing units) generally includes loading parameters for a first neuron model selected from the plurality of neuron models into a first neural processing unit, determining a first state of the first neural processing unit based at least in part on the parameters for the first neuron model, and determining a second state of the first neural processing unit based at least in part on the parameters for the first neuron model and on the first state. This method may also include updating the plurality of neuron models (e.g., by adding, deleting, or adjusting parameters for the first neuron model or another neuron model).

CLAIM OF PRIORITY UNDER 35 U.S.C. §119

This application claims benefit of U.S. Provisional Patent ApplicationSer. No. 61/728,360, filed Nov. 20, 2012 and entitled “Piecewise LinearNeuron Modeling,” U.S. Provisional Patent Application Ser. No.61/734,716, filed Dec. 7, 2012 and entitled “Piecewise Linear NeuronModeling,” U.S. Provisional Patent Application Ser. No. 61/740,633,filed Dec. 21, 2012 and entitled “Piecewise Linear Neuron Modeling,” andU.S. Provisional Patent Application Ser. No. 61/756,889, filed Jan. 25,2013 and entitled “Piecewise Linear Neuron Modeling,” all of which areherein incorporated by reference in their entireties.

BACKGROUND

1. Field

Certain aspects of the present disclosure generally relate to artificialnervous systems and, more particularly, to approximating at least aportion of a nonlinear function of a neuron model as a piecewise linearfunction and to using the resulting linearized neuron model in one ormore artificial neurons.

2. Background

An artificial neural network, which may comprise an interconnected groupof artificial neurons (i.e., neuron models), is a computational deviceor represents a method to be performed by a computational device.Artificial neural networks may have corresponding structure and/orfunction in biological neural networks. However, artificial neuralnetworks may provide innovative and useful computational techniques forcertain applications in which traditional computational techniques arecumbersome, impractical, or inadequate. Because artificial neuralnetworks can infer a function from observations, such networks areparticularly useful in applications where the complexity of the task ordata makes the design of the function by conventional techniquesburdensome.

One type of artificial neural network is the spiking neural network,which incorporates the concept of time into its operating model, as wellas neuronal and synaptic state, thereby providing a rich set ofbehaviors from which computational function can emerge in the neuralnetwork. Spiking neural networks are based on the concept that neuronsfire or “spike” at a particular time or times based on the state of theneuron, and that the time is important to neuron function. When a neuronfires, it generates a spike that travels to other neurons, which, inturn, may adjust their states based on the time this spike is received.In other words, information may be encoded in the relative or absolutetiming of spikes in the neural network.

SUMMARY

Certain aspects of the present disclosure generally relate toapproximating at least a portion of a nonlinear function of a neuronmodel as a piecewise linear function. Methods and apparatus forimplementing the resulting linearized neuron model in one or moreartificial neurons, for example, are also provided. Certain aspects ofthe present disclosure generally relate to a common and flexiblearchitecture for the implementation of the dynamics of neuron models.The design goals include low complexity, accurate modeling of thedynamics, and the ability to implement any neuron model (of one, two, ormore dimensions). The piecewise linear approximations provide a simpleway to change neuron models within such an architecture, simply bysubstituting different parameters associated with various neuron models.

Certain aspects of the present disclosure provide a method for operatingan artificial neuron. The method generally includes determining that afirst state of the artificial neuron is within a first region;determining a second state of the artificial neuron based at least inpart on a first set of linear equations, wherein the first set of linearequations is based at least in part on a first set of parameterscorresponding to the first region; determining that the second state ofthe artificial neuron is within a second region; and determining a thirdstate of the artificial neuron based at least in part on a second set oflinear equations, wherein the second set of linear equations is based atleast in part on a second set of parameters corresponding to the secondregion.

Certain aspects of the present disclosure provide an apparatus foroperating an artificial neuron. The apparatus generally includes aprocessing system and a memory coupled to the processing system. Theprocessing system is generally configured to determine that a firststate of the artificial neuron is within a first region; to determine asecond state of the artificial neuron based at least in part on a firstset of linear equations, wherein the first set of linear equations isbased at least in part on a first set of parameters corresponding to thefirst region; to determine that the second state of the artificialneuron is within a second region; and to determine a third state of theartificial neuron based at least in part on a second set of linearequations, wherein the second set of linear equations is based at leastin part on a second set of parameters corresponding to the secondregion.

Certain aspects of the present disclosure provide an apparatus foroperating an artificial neuron. The apparatus generally includes meansfor determining that a first state of the artificial neuron is within afirst region; means for determining a second state of the artificialneuron based at least in part on a first set of linear equations,wherein the first set of linear equations is based at least in part on afirst set of parameters corresponding to the first region; means fordetermining that the second state of the artificial neuron is within asecond region; and means for determining a third state of the artificialneuron based at least in part on a second set of linear equations,wherein the second set of linear equations is based at least in part ona second set of parameters corresponding to the second region.

Certain aspects of the present disclosure provide a computer programproduct for operating an artificial neuron. The computer program productgenerally includes a computer-readable medium (e.g., a storage device)having instructions executable to determine that a first state of theartificial neuron is within a first region; to determine a second stateof the artificial neuron based at least in part on a first set of linearequations, wherein the first set of linear equations is based at leastin part on a first set of parameters corresponding to the first region;to determine that the second state of the artificial neuron is within asecond region; and to determine a third state of the artificial neuronbased at least in part on a second set of linear equations, wherein thesecond set of linear equations is based at least in part on a second setof parameters corresponding to the second region.

Certain aspects of the present disclosure provide a method forimplementing a combination of a plurality of neuron models in a systemof neural processing units. The method generally includes loadingparameters for a first neuron model selected from the plurality ofneuron models into a first neural processing unit, determining a firststate of the first neural processing unit based at least in part on theparameters for the first neuron model, and determining a second state ofthe first neural processing unit based at least in part on theparameters for the first neuron model and on the first state.

Certain aspects of the present disclosure provide an apparatus forimplementing a combination of a plurality of neuron models in a systemof neural processing units. The apparatus generally includes aprocessing system and a memory coupled to the processing system. Theprocessing system is typically configured to load parameters for a firstneuron model selected from the plurality of neuron models into a firstneural processing unit, to determine a first state of the first neuralprocessing unit based at least in part on the parameters for the firstneuron model, and to determine a second state of the first neuralprocessing unit based at least in part on the parameters for the firstneuron model and on the first state.

Certain aspects of the present disclosure provide an apparatus forimplementing a combination of a plurality of neuron models in a systemof neural processing units. The apparatus generally includes means forloading parameters for a first neuron model selected from the pluralityof neuron models into a first neural processing unit, means fordetermining a first state of the first neural processing unit based atleast in part on the parameters for the first neuron model, and meansfor determining a second state of the first neural processing unit basedat least in part on the parameters for the first neuron model and on thefirst state.

Certain aspects of the present disclosure provide a computer programproduct for implementing a combination of a plurality of neuron modelsin a system of neural processing units. The computer program productgenerally includes a (non-transitory) computer-readable medium havinginstructions executable to load parameters for a first neuron modelselected from the plurality of neuron models into a first neuralprocessing unit, to determine a first state of the first neuralprocessing unit based at least in part on the parameters for the firstneuron model, and to determine a second state of the first neuralprocessing unit based at least in part on the parameters for the firstneuron model and on the first state.

Certain aspects of the present disclosure provide a method for operatingan artificial neuron. The method generally includes determining that afirst state of the artificial neuron is within a first region;determining a second state of the artificial neuron based at least inpart on a first set of linear equations, wherein the first set of linearequations is based at least in part on a first set of parameterscorresponding to the first region; determining that the second state ofthe artificial neuron is within a second region, wherein at least one ofthe first region or the second region is defined by two or moredimensions; and determining a third state of the artificial neuron basedat least in part on a second set of linear equations, wherein the secondset of linear equations is based at least in part on a second set ofparameters corresponding to the second region.

Certain aspects of the present disclosure provide an apparatus foroperating an artificial neuron. The apparatus generally includes aprocessing system and a memory coupled to the processing system. Theprocessing system is generally configured to determine that a firststate of the artificial neuron is within a first region; to determine asecond state of the artificial neuron based at least in part on a firstset of linear equations, wherein the first set of linear equations isbased at least in part on a first set of parameters corresponding to thefirst region; to determine that the second state of the artificialneuron is within a second region, wherein at least one of the firstregion or the second region is defined by two or more dimensions; and todetermine a third state of the artificial neuron based at least in parton a second set of linear equations, wherein the second set of linearequations is based at least in part on a second set of parameterscorresponding to the second region.

Certain aspects of the present disclosure provide an apparatus foroperating an artificial neuron. The apparatus generally includes meansfor determining that a first state of the artificial neuron is within afirst region; means for determining a second state of the artificialneuron based at least in part on a first set of linear equations,wherein the first set of linear equations is based at least in part on afirst set of parameters corresponding to the first region; means fordetermining that the second state of the artificial neuron is within asecond region, wherein at least one of the first region or the secondregion is defined by two or more dimensions; and means for determining athird state of the artificial neuron based at least in part on a secondset of linear equations, wherein the second set of linear equations isbased at least in part on a second set of parameters corresponding tothe second region.

Certain aspects of the present disclosure provide a computer programproduct for operating an artificial neuron. The computer program productgenerally includes a computer-readable medium (e.g., a storage device orother non-transitory medium) having instructions executable to determinethat a first state of the artificial neuron is within a first region; todetermine a second state of the artificial neuron based at least in parton a first set of linear equations, wherein the first set of linearequations is based at least in part on a first set of parameterscorresponding to the first region; to determine that the second state ofthe artificial neuron is within a second region, wherein at least one ofthe first region or the second region is defined by two or moredimensions; and to determine a third state of the artificial neuronbased at least in part on a second set of linear equations, wherein thesecond set of linear equations is based at least in part on a second setof parameters corresponding to the second region.

BRIEF DESCRIPTION OF THE DRAWINGS

So that the manner in which the above-recited features of the presentdisclosure can be understood in detail, a more particular description,briefly summarized above, may be had by reference to aspects, some ofwhich are illustrated in the appended drawings. It is to be noted,however, that the appended drawings illustrate only certain typicalaspects of this disclosure and are therefore not to be consideredlimiting of its scope, for the description may admit to other equallyeffective aspects.

FIG. 1 illustrates an example network of neurons in accordance withcertain aspects of the present disclosure.

FIG. 2 illustrates an example processing unit (neuron) of acomputational network (neural system or neural network), in accordancewith certain aspects of the present disclosure.

FIG. 3 illustrates an example spike-timing dependent plasticity (STDP)curve in accordance with certain aspects of the present disclosure.

FIG. 4 is an example graph of state for an artificial neuron,illustrating a positive regime and a negative regime for definingbehavior of the neuron, in accordance with certain aspects of thepresent disclosure.

FIGS. 5A and 5B illustrate example graphs of membrane voltage v andrecovery current u, respectively, versus time for comparing thenonlinear time-varying simple model to an example linearization based ona Taylor expansion method, in accordance with certain aspects of thepresent disclosure.

FIGS. 6A and 6B illustrate example graphs of membrane voltage v andrecovery current u, respectively, versus time for the subthresholddynamics of the Hunzinger Cold model, in accordance with certain aspectsof the present disclosure.

FIG. 7 illustrates a summary of various approaches to piecewise linearneuron modeling, in accordance with certain aspects of the presentdisclosure.

FIG. 8 illustrates an example of piecewise linear neuron modeling withpartitioning in terms of a single dimension, in accordance with certainaspects of the present disclosure.

FIG. 9 illustrates an example of generalized linear neuron modeling withthree rectangular regions defined by two dimensions, in accordance withcertain aspects of the present disclosure.

FIG. 10 illustrates an example of generalized linear neuron modelingwith four regions defined by two dimensions, which may be used toexhibit decaying, sustained, and growing subthreshold oscillationbehavior, in accordance with certain aspects of the present disclosure.

FIG. 11 illustrates an example of generalized linear neuron modelingwith five regions having varying shapes, in accordance with certainaspects of the present disclosure.

FIGS. 12A-C illustrate example plots of decaying, sustained, and growingsubthreshold oscillations, respectively, in accordance with certainaspects of the present disclosure.

FIG. 13 illustrates an example of generalized linear neuron modelingwith six regions defined by two dimensions, which may be used formulti-stage decay, multi-stage sustained, and/or multi-stage growthregions to support gradual decaying, multiple sustained oscillations,and/or gradual growing of subthreshold oscillations, in accordance withcertain aspects of the present disclosure.

FIG. 14 is a flow diagram of example operations for operating anartificial neuron, in accordance with certain aspects of the presentdisclosure.

FIG. 14A illustrates example means capable of performing the operationsshown in FIG. 14.

FIGS. 15A-D illustrate implementation of a common and flexible neuralarchitecture for a single neural processing unit, where parameters forneuron models can be selected, loaded, accessed, added, deleted, and/orupdated, in accordance with certain aspects of the present disclosure.

FIG. 16 is a flow diagram of example operations for implementing acombination of a plurality of neuron models in a system of neuralprocessing units, in accordance with certain aspects of the presentdisclosure.

FIG. 16A illustrates example means capable of performing the operationsshown in FIG. 16.

FIG. 17 illustrates an example implementation for determining states ofan artificial neuron using a general-purpose processor, in accordancewith certain aspects of the present disclosure.

FIG. 18 illustrates an example implementation for determining states ofan artificial neuron where a memory may be interfaced with individualdistributed processing units, in accordance with certain aspects of thepresent disclosure.

FIG. 19 illustrates an example implementation for determining states ofan artificial neuron based on distributed memories and distributedprocessing units, in accordance with certain aspects of the presentdisclosure.

FIG. 20 illustrates an example implementation of a neural network inaccordance with certain aspects of the present disclosure.

FIG. 21 is a block diagram of an example implementation of piecewiselinear neuron modeling in which parameters used to update states of anartificial neuron are fetched from memory according to the quantizationof the current state, in accordance with certain aspects of the presentdisclosure.

DETAILED DESCRIPTION

Various aspects of the disclosure are described more fully hereinafterwith reference to the accompanying drawings. This disclosure may,however, be embodied in many different forms and should not be construedas limited to any specific structure or function presented throughoutthis disclosure. Rather, these aspects are provided so that thisdisclosure will be thorough and complete, and will fully convey thescope of the disclosure to those skilled in the art. Based on theteachings herein one skilled in the art should appreciate that the scopeof the disclosure is intended to cover any aspect of the disclosuredisclosed herein, whether implemented independently of or combined withany other aspect of the disclosure. For example, an apparatus may beimplemented or a method may be practiced using any number of the aspectsset forth herein. In addition, the scope of the disclosure is intendedto cover such an apparatus or method which is practiced using otherstructure, functionality, or structure and functionality in addition toor other than the various aspects of the disclosure set forth herein. Itshould be understood that any aspect of the disclosure disclosed hereinmay be embodied by one or more elements of a claim.

The word “exemplary” is used herein to mean “serving as an example,instance, or illustration.” Any aspect described herein as “exemplary”is not necessarily to be construed as preferred or advantageous overother aspects.

Although particular aspects are described herein, many variations andpermutations of these aspects fall within the scope of the disclosure.Although some benefits and advantages of the preferred aspects arementioned, the scope of the disclosure is not intended to be limited toparticular benefits, uses or objectives. Rather, aspects of thedisclosure are intended to be broadly applicable to differenttechnologies, system configurations, networks and protocols, some ofwhich are illustrated by way of example in the figures and in thefollowing description of the preferred aspects. The detailed descriptionand drawings are merely illustrative of the disclosure rather thanlimiting, the scope of the disclosure being defined by the appendedclaims and equivalents thereof.

An Example Neural System

FIG. 1 illustrates an example neural system 100 with multiple levels ofneurons in accordance with certain aspects of the present disclosure.The neural system 100 may comprise a level of neurons 102 connected toanother level of neurons 106 though a network of synaptic connections104 (i.e., feed-forward connections). For simplicity, only two levels ofneurons are illustrated in FIG. 1, although fewer or more levels ofneurons may exist in a typical neural system. It should be noted thatsome of the neurons may connect to other neurons of the same layerthrough lateral connections. Furthermore, some of the neurons mayconnect back to a neuron of a previous layer through feedbackconnections.

As illustrated in FIG. 1, each neuron in the level 102 may receive aninput signal 108 that may be generated by a plurality of neurons of aprevious level (not shown in FIG. 1). The signal 108 may represent aninput (e.g., an input current) to the level 102 neuron. Such inputs maybe accumulated on the neuron membrane to charge a membrane potential.When the membrane potential reaches its threshold value, the neuron mayfire and generate an output spike to be transferred to the next level ofneurons (e.g., the level 106). Such behavior can be emulated orsimulated in hardware and/or software, including analog and digitalimplementations.

In biological neurons, the output spike generated when a neuron fires isreferred to as an action potential. This electrical signal is arelatively rapid, transient, all-or nothing nerve impulse, having anamplitude of roughly 100 mV and a duration of about 1 ms. In aparticular aspect of a neural system having a series of connectedneurons (e.g., the transfer of spikes from one level of neurons toanother in FIG. 1), every action potential has basically the sameamplitude and duration, and thus, the information in the signal isrepresented only by the frequency and number of spikes (or the time ofspikes), not by the amplitude. The information carried by an actionpotential is determined by the spike, the neuron that spiked, and thetime of the spike relative to one or more other spikes.

The transfer of spikes from one level of neurons to another may beachieved through the network of synaptic connections (or simply“synapses”) 104, as illustrated in FIG. 1. The synapses 104 may receiveoutput signals (i.e., spikes) from the level 102 neurons (pre-synapticneurons relative to the synapses 104). For certain aspects, thesesignals may be scaled according to adjustable synaptic weights w₁^((i,i+1)), . . . , w_(P) ^((i,i+1)) (where P is a total number ofsynaptic connections between the neurons of levels 102 and 106). Forother aspects, the synapses 104 may not apply any synaptic weights.Further, the (scaled) signals may be combined as an input signal of eachneuron in the level 106 (post-synaptic neurons relative to the synapses104). Every neuron in the level 106 may generate output spikes 110 basedon the corresponding combined input signal. The output spikes 110 may bethen transferred to another level of neurons using another network ofsynaptic connections (not shown in FIG. 1).

Biological synapses may be classified as either electrical or chemical.While electrical synapses are used primarily to send excitatory signals,chemical synapses can mediate either excitatory or inhibitory(hyperpolarizing) actions in postsynaptic neurons and can also serve toamplify neuronal signals. Excitatory signals typically depolarize themembrane potential (i.e., increase the membrane potential with respectto the resting potential). If enough excitatory signals are receivedwithin a certain period to depolarize the membrane potential above athreshold, an action potential occurs in the postsynaptic neuron. Incontrast, inhibitory signals generally hyperpolarize (i.e., lower) themembrane potential Inhibitory signals, if strong enough, can counteractthe sum of excitatory signals and prevent the membrane potential fromreaching threshold. In addition to counteracting synaptic excitation,synaptic inhibition can exert powerful control over spontaneously activeneurons. A spontaneously active neuron refers to a neuron that spikeswithout further input, for example, due to its dynamics or feedback. Bysuppressing the spontaneous generation of action potentials in theseneurons, synaptic inhibition can shape the pattern of firing in aneuron, which is generally referred to as sculpturing. The varioussynapses 104 may act as any combination of excitatory or inhibitorysynapses, depending on the behavior desired.

The neural system 100 may be emulated by a general purpose processor, adigital signal processor (DSP), an application specific integratedcircuit (ASIC), a field programmable gate array (FPGA) or otherprogrammable logic device (PLD), discrete gate or transistor logic,discrete hardware components, a software module executed by a processor,or any combination thereof. The neural system 100 may be utilized in alarge range of applications, such as image and pattern recognition,machine learning, motor control, and the like. Each neuron (or neuronmodel) in the neural system 100 may be implemented as a neuron circuit.The neuron membrane charged to the threshold value initiating the outputspike may be implemented, for example, as a capacitor that integrates anelectrical current flowing through it.

In an aspect, the capacitor may be eliminated as the electrical currentintegrating device of the neuron circuit, and a smaller memristorelement may be used in its place. This approach may be applied in neuroncircuits, as well as in various other applications where bulkycapacitors are utilized as electrical current integrators. In addition,each of the synapses 104 may be implemented based on a memristorelement, wherein synaptic weight changes may relate to changes of thememristor resistance. With nanometer feature-sized memristors, the areaof neuron circuit and synapses may be substantially reduced, which maymake implementation of a very large-scale neural system hardwareimplementation practical.

Functionality of a neural processor that emulates the neural system 100may depend on weights of synaptic connections, which may controlstrengths of connections between neurons. The synaptic weights may bestored in a non-volatile memory in order to preserve functionality ofthe processor after being powered down. In an aspect, the synapticweight memory may be implemented on a separate external chip from themain neural processor chip. The synaptic weight memory may be packagedseparately from the neural processor chip as a replaceable memory card.This may provide diverse functionalities to the neural processor,wherein a particular functionality may be based on synaptic weightsstored in a memory card currently attached to the neural processor.

FIG. 2 illustrates an example 200 of a processing unit (e.g., anartificial neuron 202) of a computational network (e.g., a neural systemor a neural network) in accordance with certain aspects of the presentdisclosure. For example, the neuron 202 may correspond to any of theneurons of levels 102 and 106 from FIG. 1. The neuron 202 may receivemultiple input signals 204 ₁-204 _(N) (x₁-x_(N)), which may be signalsexternal to the neural system, or signals generated by other neurons ofthe same neural system, or both. The input signal may be a current or avoltage, real-valued or complex-valued. The input signal may comprise anumerical value with a fixed-point or a floating-point representation.These input signals may be delivered to the neuron 202 through synapticconnections that scale the signals according to adjustable synapticweights 206 ₁-206 _(N) (w₁-w_(N)), where N may be a total number ofinput connections of the neuron 202.

The neuron 202 may combine the scaled input signals and use the combinedscaled inputs to generate an output signal 208 (i.e., a signal y). Theoutput signal 208 may be a current, or a voltage, real-valued orcomplex-valued. The output signal may comprise a numerical value with afixed-point or a floating-point representation. The output signal 208may be then transferred as an input signal to other neurons of the sameneural system, or as an input signal to the same neuron 202, or as anoutput of the neural system.

The processing unit (neuron) 202 may be emulated by an electricalcircuit, and its input and output connections may be emulated by wireswith synaptic circuits. The processing unit 202, its input and outputconnections may also be emulated by a software code. The processing unit202 may also be emulated by an electric circuit, whereas its input andoutput connections may be emulated by a software code. In an aspect, theprocessing unit 202 in the computational network may comprise an analogelectrical circuit. In another aspect, the processing unit 202 maycomprise a digital electrical circuit. In yet another aspect, theprocessing unit 202 may comprise a mixed-signal electrical circuit withboth analog and digital components. The computational network maycomprise processing units in any of the aforementioned forms. Thecomputational network (neural system or neural network) using suchprocessing units may be utilized in a large range of applications, suchas image and pattern recognition, machine learning, motor control, andthe like.

During the course of training a neural network, synaptic weights (e.g.,the weights w₁ ^((i,i+1)), . . . , w_(P) ^((i,i+1)) from FIG. 1 and/orthe weights 206 ₁-206 _(N) from FIG. 2) may be initialized with randomvalues and increased or decreased according to a learning rule. Someexamples of the learning rule are the spike-timing-dependent plasticity(STDP) learning rule, the Hebb rule, the Oja rule, theBienenstock-Copper-Munro (BCM) rule, etc. Very often, the weights maysettle to one of two values (i.e., a bimodal distribution of weights).This effect can be utilized to reduce the number of bits per synapticweight, increase the speed of reading and writing from/to a memorystoring the synaptic weights, and to reduce power consumption of thesynaptic memory.

Synapse Type

In hardware and software models of neural networks, processing ofsynapse related functions can be based on synaptic type. Synapse typesmay comprise non-plastic synapses (no changes of weight and delay),plastic synapses (weight may change), structural delay plastic synapses(weight and delay may change), fully plastic synapses (weight, delay andconnectivity may change), and variations thereupon (e.g., delay maychange, but no change in weight or connectivity). The advantage of thisis that processing can be subdivided. For example, non-plastic synapsesmay not require plasticity functions to be executed (or waiting for suchfunctions to complete). Similarly, delay and weight plasticity may besubdivided into operations that may operate in together or separately,in sequence or in parallel. Different types of synapses may havedifferent lookup tables or formulas and parameters for each of thedifferent plasticity types that apply. Thus, the methods would accessthe relevant tables for the synapse's type.

There are further implications of the fact that spike-timing dependentstructural plasticity may be executed independently of synapticplasticity. Structural plasticity may be executed even if there is nochange to weight magnitude (e.g., if the weight has reached a minimum ormaximum value, or it is not changed due to some other reason) sincestructural plasticity (i.e., an amount of delay change) may be a directfunction of pre-post spike time difference. Alternatively, it may be setas a function of the weight change amount or based on conditionsrelating to bounds of the weights or weight changes. For example, asynaptic delay may change only when a weight change occurs or if weightsreach zero, but not if the weights are maxed out. However, it can beadvantageous to have independent functions so that these processes canbe parallelized reducing the number and overlap of memory accesses.

Determination of Synaptic Plasticity

Neuroplasticity (or simply “plasticity”) is the capacity of neurons andneural networks in the brain to change their synaptic connections andbehavior in response to new information, sensory stimulation,development, damage, or dysfunction. Plasticity is important to learningand memory in biology, as well as to computational neuroscience andneural networks. Various forms of plasticity have been studied, such assynaptic plasticity (e.g., according to the Hebbian theory),spike-timing-dependent plasticity (STDP), non-synaptic plasticity,activity-dependent plasticity, structural plasticity, and homeostaticplasticity.

STDP is a learning process that adjusts the strength of synapticconnections between neurons, such as those in the brain. The connectionstrengths are adjusted based on the relative timing of a particularneuron's output and received input spikes (i.e., action potentials).Under the STDP process, long-term potentiation (LTP) may occur if aninput spike to a certain neuron tends, on average, to occur immediatelybefore that neuron's output spike. Then, that particular input is madesomewhat stronger. In contrast, long-term depression (LTD) may occur ifan input spike tends, on average, to occur immediately after an outputspike. Then, that particular input is made somewhat weaker, hence thename “spike-timing-dependent plasticity.” Consequently, inputs thatmight be the cause of the post-synaptic neuron's excitation are madeeven more likely to contribute in the future, whereas inputs that arenot the cause of the post-synaptic spike are made less likely tocontribute in the future. The process continues until a subset of theinitial set of connections remains, while the influence of all others isreduced to zero or near zero.

Since a neuron generally produces an output spike when many of itsinputs occur within a brief period (i.e., being sufficiently cumulativeto cause the output,), the subset of inputs that typically remainsincludes those that tended to be correlated in time. In addition, sincethe inputs that occur before the output spike are strengthened, theinputs that provide the earliest sufficiently cumulative indication ofcorrelation will eventually become the final input to the neuron.

The STDP learning rule may effectively adapt a synaptic weight of asynapse connecting a pre-synaptic neuron to a post-synaptic neuron as afunction of time difference between spike time t_(pre) of thepre-synaptic neuron and spike time t_(post) of the post-synaptic neuron(i.e., t=t_(post)−t_(pre)). A typical formulation of the STDP is toincrease the synaptic weight (i.e., potentiate the synapse) if the timedifference is positive (the pre-synaptic neuron fires before thepost-synaptic neuron), and decrease the synaptic weight (i.e., depressthe synapse) if the time difference is negative (the post-synapticneuron fires before the pre-synaptic neuron).

In the STDP process, a change of the synaptic weight over time may betypically achieved using an exponential decay, as given by,

$\begin{matrix}{{\Delta \; {w(t)}} = \left\{ \begin{matrix}{{{a_{+}^{{- t}/k_{+}}} + \mu},{t > 0}} \\{{a_{-}^{t/k_{-}}},{t < 0},}\end{matrix} \right.} & (1)\end{matrix}$

where k₊ and k⁻ are time constants for positive and negative timedifference, respectively, a₊ and a⁻ are corresponding scalingmagnitudes, and μ is an offset that may be applied to the positive timedifference and/or the negative time difference.

FIG. 3 illustrates an example graph diagram 300 of a synaptic weightchange as a function of relative timing of pre-synaptic andpost-synaptic spikes in accordance with STDP. If a pre-synaptic neuronfires before a post-synaptic neuron, then a corresponding synapticweight may be increased, as illustrated in a portion 302 of the graph300. This weight increase can be referred to as an LTP of the synapse.It can be observed from the graph portion 302 that the amount of LTP maydecrease roughly exponentially as a function of the difference betweenpre-synaptic and post-synaptic spike times. The reverse order of firingmay reduce the synaptic weight, as illustrated in a portion 304 of thegraph 300, causing an LTD of the synapse.

As illustrated in the graph 300 in FIG. 3, a negative offset may beapplied to the LTP (causal) portion 302 of the STDP graph. A point ofcross-over 306 of the x-axis (y=0) may be configured to coincide withthe maximum time lag for considering correlation for causal inputs fromlayer i−1. In the case of a frame-based input (i.e., an input is in theform of a frame of a particular duration comprising spikes or pulses),the offset value μ can be computed to reflect the frame boundary. Afirst input spike (pulse) in the frame may be considered to decay overtime either as modeled by a post-synaptic potential directly or in termsof the effect on neural state. If a second input spike (pulse) in theframe is considered correlated or relevant of a particular time frame,then the relevant times before and after the frame may be separated atthat time frame boundary and treated differently in plasticity terms byoffsetting one or more parts of the STDP curve such that the value inthe relevant times may be different (e.g., negative for greater than oneframe and positive for less than one frame). For example, the negativeoffset μ may be set to offset LTP such that the curve actually goesbelow zero at a pre-post time greater than the frame time and it is thuspart of LTD instead of LTP.

Neuron Models and Operation

There are some general principles for designing a useful spiking neuronmodel. A good neuron model may have rich potential behavior in terms oftwo computational regimes: coincidence detection and functionalcomputation. Moreover, a good neuron model should have two elements toallow temporal coding: arrival time of inputs affects output time andcoincidence detection can have a narrow time window. Finally, to becomputationally attractive, a good neuron model may have a closed-formsolution in continuous time and have stable behavior including nearattractors and saddle points. In other words, a useful neuron model isone that is practical and that can be used to model rich, realistic andbiologically-consistent behaviors, as well as be used to both engineerand reverse engineer neural circuits.

A neuron model may depend on events, such as an input arrival, outputspike or other event whether internal or external. To achieve a richbehavioral repertoire, a state machine that can exhibit complexbehaviors may be desired. If the occurrence of an event itself, separatefrom the input contribution (if any) can influence the state machine andconstrain dynamics subsequent to the event, then the future state of thesystem is not only a function of a state and input, but rather afunction of a state, event, and input.

In an aspect, a neuron n may be modeled as a spikingleaky-integrate-and-fire neuron with a membrane voltage v_(n)(t)governed by the following dynamics,

$\begin{matrix}{{\frac{{v_{n}(t)}}{t} = {{\alpha \; {v_{n}(t)}} + {\beta {\sum\limits_{m}{w_{m,n}{y_{m}\left( {t - {\Delta \; t_{m,n}}} \right)}}}}}},} & (2)\end{matrix}$

where α and β are parameters, w_(m,n) is a synaptic weight for thesynapse connecting a pre-synaptic neuron m to a post-synaptic neuron n,and y_(m)(t) is the spiking output of the neuron m that may be delayedby dendritic or axonal delay according to Δt_(m,n) until arrival at theneuron n's soma.

It should be noted that there is a delay from the time when sufficientinput to a post-synaptic neuron is established until the time when thepost-synaptic neuron actually fires. In a dynamic spiking neuron model,such as Izhikevich's simple model, a time delay may be incurred if thereis a difference between a depolarization threshold v_(t) and a peakspike voltage v_(peak). For example, in the simple model, neuron somadynamics can be governed by the pair of differential equations forvoltage and recovery, i.e.,

$\begin{matrix}{{\frac{v}{t} = {\left( {{{k\left( {v - v_{t}} \right)}\left( {v - v_{r}} \right)} - u + I} \right)/C}},} & (3) \\{\frac{u}{t} = {{a\left( {{b\left( {v - v_{r}} \right)} - u} \right)}.}} & (4)\end{matrix}$

where v is a membrane potential, u is a membrane recovery variable, k isa parameter that describes time scale of the membrane potential v, a isa parameter that describes time scale of the recovery variable u, b is aparameter that describes sensitivity of the recovery variable u to thesubthreshold fluctuations of the membrane potential v, v_(r) is amembrane resting potential, I is a synaptic current, and C is amembrane's capacitance. In accordance with this model, the neuron isdefined to spike when v>v_(peak).

Hunzinger Cold Model

The Hunzinger Cold neuron model is a minimal dual-regime spiking lineardynamical model that can reproduce a rich variety of neural behaviors.The model's one- or two-dimensional linear dynamics can have tworegimes, wherein the time constant (and coupling) can depend on theregime. In the subthreshold regime, the time constant, negative byconvention, represents leaky channel dynamics generally acting to returna cell to rest in biologically-consistent linear fashion. The timeconstant in the supra-threshold regime, positive by convention, reflectsanti-leaky channel dynamics generally driving a cell to spike whileincurring latency in spike-generation.

As illustrated in FIG. 4, the dynamics of the model may be divided intotwo (or more) regimes. These regimes may be called the negative regime402 (also interchangeably referred to as the leaky-integrate-and-fire(LIF) regime, not to be confused with the LIF neuron model) and thepositive regime 404 (also interchangeably referred to as theanti-leaky-integrate-and-fire (ALIF) regime, not to be confused with theALIF neuron model). In the negative regime 402, the state tends towardrest (v⁻) at the time of a future event. In this negative regime, themodel generally exhibits temporal input detection properties and othersubthreshold behavior. In the positive regime 404, the state tendstoward a spiking event (v_(s)). In this positive regime, the modelexhibits computational properties, such as incurring a latency to spikedepending on subsequent input events. Formulation of dynamics in termsof events and separation of the dynamics into these two regimes arefundamental characteristics of the model.

Linear dual-regime bi-dimensional dynamics (for states v and u) may bedefined by convention as,

$\begin{matrix}{{\tau_{\rho}\frac{v}{t}} = {v + q_{\rho}}} & (5) \\{{{- \tau_{u}}\frac{u}{t}} = {u + r}} & (6)\end{matrix}$

where q_(ρ) and r are the linear transformation variables for coupling.

The symbol ρ is used herein to denote the dynamics regime with theconvention to replace the symbol ρ with the sign “−” or “+” for thenegative and positive regimes, respectively, when discussing orexpressing a relation for a specific regime.

The model state is defined by a membrane potential (voltage) v andrecovery current u. In basic form, the regime is essentially determinedby the model state. There are subtle, but important aspects of theprecise and general definition, but for the moment, consider the modelto be in the positive regime 404 if the voltage v is above a threshold(v₊) and otherwise in the negative regime 402.

The regime-dependent time constants include τ⁻ which is the negativeregime time constant, and τ₊ which is the positive regime time constant.The recovery current time constant is typically independent of regime.For convenience, the negative regime time constant τ⁻ is typicallyspecified as a negative quantity to reflect decay so that the sameexpression for voltage evolution may be used as for the positive regimein which the exponent and τ₊ will generally be positive, as will beτ_(u).

The dynamics of the two state elements may be coupled at events bytransformations offsetting the states from their null-clines, where thetransformation variables are

q _(ρ)=−τ_(ρ) /βu−v _(ρ)  (7)

r=δ(v+ε)  (8)

where δ, ε, β and v⁻, v₊ are parameters. The two values for v_(ρ) arethe base for reference voltages for the two regimes. The parameter v⁻ isthe base voltage for the negative regime, and the membrane potentialwill generally decay toward v⁻ in the negative regime. The parameter v₊is the base voltage for the positive regime, and the membrane potentialwill generally tend away from v₊ in the positive regime.

The null-clines for v and u are given by the negative of thetransformation variables q_(ρ) and r, respectively. The parameter δ is ascale factor controlling the slope of the u null-cline. The parameter εis typically set equal to −v⁻. The parameter β is a resistance valuecontrolling the slope of the v null-clines in both regimes. The τ_(ρ)time-constant parameters control not only the exponential decays, butalso the null-cline slopes in each regime separately.

The model is defined to spike when the voltage v reaches a value v_(s).Subsequently, the state is typically reset at a reset event (whichtechnically may be one and the same as the spike event):

v={circumflex over (v)} ⁻  (9)

u=u+Δu  (10)

where {circumflex over (v)}⁻ and Δu are parameters. The reset voltage{circumflex over (v)}⁻ is typically set to v⁻.

By a principle of momentary coupling, a closed form solution is possiblenot only for state (and with a single exponential term), but also forthe time required to reach a particular state. The close form statesolutions are

$\begin{matrix}{{v\left( {t + {\Delta \; t}} \right)} = {{\left( {{v(t)} + q_{\rho}} \right)^{\frac{\Delta \; t}{\tau_{\rho}}}} - q_{\rho}}} & (11) \\{{u\left( {t + {\Delta \; t}} \right)} = {{\left( {{u(t)} + r} \right)^{- \frac{\Delta \; t}{\tau_{u}}}} - r}} & (12)\end{matrix}$

Therefore, the model state may be updated only upon events such as uponan input (pre-synaptic spike) or output (post-synaptic spike).Operations may also be performed at any particular time (whether or notthere is input or output).

Moreover, by the momentary coupling principle, the time of apost-synaptic spike may be anticipated so the time to reach a particularstate may be determined in advance without iterative techniques orNumerical Methods (e.g., the Euler numerical method). Given a priorvoltage state v₀, the time delay until voltage state v_(f) is reached isgiven by

$\begin{matrix}{{\Delta \; t} = {\tau_{\rho}\log \; \frac{v_{f} + q_{\rho}}{v_{0} + q_{\rho}}}} & (13)\end{matrix}$

If a spike is defined as occurring at the time the voltage state vreaches v_(s), then the closed-form solution for the amount of time, orrelative delay, until a spike occurs as measured from the time that thevoltage is at a given state v is

$\begin{matrix}{{\Delta \; t_{S}} = \left\{ \begin{matrix}{\tau_{+}\log \; \frac{v_{S} + q_{+}}{v + q_{+}}} & {{{if}\mspace{14mu} v} > {\hat{v}}_{+}} \\\infty & {otherwise}\end{matrix} \right.} & (14)\end{matrix}$

where {circumflex over (v)}₊ is typically set to parameter v₊, althoughother variations may be possible.

The above definitions of the model dynamics depend on whether the modelis in the positive or negative regime. As mentioned, the coupling andthe regime ρ may be computed upon events. For purposes of statepropagation, the regime and coupling (transformation) variables may bedefined based on the state at the time of the last (prior) event. Forpurposes of subsequently anticipating spike output time, the regime andcoupling variable may be defined based on the state at the time of thenext (current) event.

There are several possible implementations of the Cold model, andexecuting the simulation, emulation or model in time. This includes, forexample, event-update, step-event update, and step-update modes. Anevent update is an update where states are updated based on events or“event update” (at particular moments). A step update is an update whenthe model is updated at intervals (e.g., 1 ms). This does notnecessarily require iterative methods or Numerical methods. Anevent-based implementation is also possible at a limited time resolutionin a step-based simulator by only updating the model if an event occursat or between steps or by “step-event” update.

Neural Coding

A useful neural network model, such as one composed of the artificialneurons 102, 106 of FIG. 1, may encode information via any of varioussuitable neural coding schemes, such as coincidence coding, temporalcoding or rate coding. In coincidence coding, information is encoded inthe coincidence (or temporal proximity) of action potentials (spikingactivity) of a neuron population. In temporal coding, a neuron encodesinformation through the precise timing of action potentials (i.e.,spikes) whether in absolute time or relative time. Information may thusbe encoded in the relative timing of spikes among a population ofneurons. In contrast, rate coding involves coding the neural informationin the firing rate or population firing rate.

If a neuron model can perform temporal coding, then it can also performrate coding (since rate is just a function of timing or inter-spikeintervals). To provide for temporal coding, a good neuron model shouldhave two elements: (1) arrival time of inputs affects output time; and(2) coincidence detection can have a narrow time window. Connectiondelays provide one means to expand coincidence detection to temporalpattern decoding because by appropriately delaying elements of atemporal pattern, the elements may be brought into timing coincidence.

Arrival Time

In a good neuron model, the time of arrival of an input should have aneffect on the time of output. A synaptic input—whether a Dirac deltafunction or a shaped post-synaptic potential (PSP), whether excitatory(EPSP) or inhibitory (IPSP)—has a time of arrival (e.g., the time of thedelta function or the start or peak of a step or other input function),which may be referred to as the input time. A neuron output (i.e., aspike) has a time of occurrence (wherever it is measured, e.g., at thesoma, at a point along the axon, or at an end of the axon), which may bereferred to as the output time. That output time may be the time of thepeak of the spike, the start of the spike, or any other time in relationto the output waveform. The overarching principle is that the outputtime depends on the input time.

One might at first glance think that all neuron models conform to thisprinciple, but this is generally not true. For example, rate-basedmodels do not have this feature. Many spiking models also do notgenerally conform. A leaky-integrate-and-fire (LIF) model does not fireany faster if there are extra inputs (beyond threshold). Moreover,models that might conform if modeled at very high timing resolutionoften will not conform when timing resolution is limited, such as to 1ms steps.

Inputs

An input to a neuron model may include Dirac delta functions, such asinputs as currents, or conductance-based inputs. In the latter case, thecontribution to a neuron state may be continuous or state-dependent.

Example Piecewise Linear Neuron Modeling

Mathematical models for the dynamics of a neuron have been sought andstudied for decades. A variety of neuron models have been proposed,varying in complexity and the accuracy with which the models match theirbiological counterparts. Fundamentally, all neuron models attempt tocapture the nonlinear behavior of a cell membrane voltage due to theinteractions of a large variety of ion channels and have a commonstarting point, namely the mathematical description provided by thebreakthrough work of Hodgkin-Huxley in the 1950s.

Over the years, neuroscientists have converged, in large part, towardstwo-dimensional neuron models, which appear to provide a good tradeoffbetween the ability of duplicating the measured behavior of thebiological cells they seek to model and the ease and speed with whichthey can be analyzed and simulated. The most common two-dimensionalmodels, all formulated via a pair of differential equations, aredescribed below.

However, conventional neuron model implementations lack flexibility.Typically, implementation of a neuron model for an artificial neuron,for example, is predicated on selecting a particular neuron modelbeforehand. The downside to this approach is the difficulty ofimplementing a completely different or even a slightly modified neuronmodel.

Certain aspects of the present disclosure take advantage of the factthat the differential equations for various dynamical neuron models maybe considered as equivalent, but for the nonlinear function F(v) tryingto capture the nonlinear behavior of the cell membrane voltage, asdescribed above. With this realization, one approach includeslinearizing, at finite quantization intervals, the nonlinear functionsin the differential equations used to model the dynamics of the neuronsto provide a piecewise linear approximation. Advantages of such anapproach include the ability to derive solutions to the dynamics ineither continuous time or discrete time with relative ease and a generalmathematical framework whereby any neuron model can be analyzed andsimulated. These solutions provide parameters corresponding to eachquantization interval for a given neuron model, such that implementationof a different neuron model for certain aspects of the presentdisclosure may involve a simple substitution of the parameters.

Certain aspects of the present disclosure apply this piecewiselinearization approach to a functional that includes the synapticcurrent. This more generalized approach leads to a system matrix whichis a function, among other variables, of the time-varying synapticconductance. The present disclosure first investigates some approximatesolutions to the resulting piecewise linear time-varying system. Next,the present disclosure investigates a piecewise linear time-invariant(LTI) system obtained from the previous time-varying system byapproximating the time-varying conductance by a constant over a giventime interval. This more precise approach, with a functional thatincludes the synaptic current, leads to system matrices with a largedynamic range to be covered by the quantization process and therefore torelatively large memory demands to store various pre-calculatedmatrices. The benefit of this approach is a more accurate approximationof the true dynamics of the original nonlinear time-varying system. Anapproach to alleviate the memory demands by computing the relevantmatrices, rather than storing them, over a portion of the dynamic rangeis also presented.

Certain aspects of the present disclosure provide a means to realize acommon architecture that supports any one-dimensional, two-dimensional,or higher dimensioned neuron models. With this flexible architecture,any of various suitable neuron models may be executed and substituted asdesired. For example, the neuron model may include at least one of anIzhikevich simple model, an exponential-integrate-and-fire (EIF) model,a FitzHugh-Nagumo model, a quartic model, or a Hunzinger Cold model asdescribed above and in U.S. patent application Ser. No. 13/483,811[Atty. Dkt. No. 122024], entitled “Dynamical Event Neuron and SynapseModels for Learning Spiking Neural Networks” and filed May 30, 2012,herein incorporated by reference. Such neuron models may be implementedusing the piecewise linear approximation described herein.

INTRODUCTION

One place to begin is with the description of the differential equationsmodeling neuron dynamics. Although the description focuses ontwo-dimensional neuron models, the approach can be extended to higherdimensional models or applied to one-dimensional models, as well.

$\begin{matrix}{{C\frac{}{t}v} = {{F(v)} - u + I}} & (15) \\{{\frac{}{t}u} = {a\left( {{b \cdot \left( {v - v_{r}} \right)} - u} \right)}} & (16)\end{matrix}$

The above equations qualitatively describe the dynamics of an artificialneuron (for notational simplicity, the time dependency in the variablesis omitted). These equations are the result of simplifications of theHodgkin-Huxley four dimensional model to two dimensions, represented bythe variables v and u. The variable v captures the behavior of theneuron membrane voltage and sodium activation, while u represents an“accommodation” or “recovery” variable attempting to capture the slowerbehavior of the potassium activation and sodium inactivation, thusreducing the four variables of the Hodgkin-Huxley model to two. Thevariable I in Eq. (15) represents the input current. A more generictwo-dimensional model may be of the form

${C\frac{}{t}v} = {{F\left( {v,u} \right)} + I}$${\frac{}{t}u} = {G\left( {v,u} \right)}$

where both differential equations may contain nonlinear terms. Here, thefocus is on models such as those described by Eqs. (15) and (16), butthe same methodology developed throughout can be applied to the abovedescription, as well.

Fundamentally, the most popular neuron models that have been suggestedin the literature differ by the choice of the function F(v) in Eq. (15).Some examples include a quadratic function, as suggested by Izhikevich(also referred as the simple model):

F(v)=k(v−v _(r))(v−v _(t))  (17)

a linear-plus-exponential function, per Brette and Gerstner:

$\begin{matrix}{{F(v)} = {{g_{L} \cdot \Delta \cdot ^{\frac{v - v_{t}}{\Delta}}} + {g_{L} \cdot \left( {E_{L} - v} \right)}}} & (18)\end{matrix}$

a linear-plus-cubic function representing the FitzHugh-Nagumo model:

$\begin{matrix}{{F(v)} = {v - {\frac{1}{3}v^{3}}}} & (19)\end{matrix}$

a linear-plus-quartic term, referred to as the quartic model, perTouboul and Gerstner

F(v)=2a·v+v ⁴  (20)

and finally what may be referred to as “intrinsic conductance” models,defined as

F(v)=G(v)·v+p(v)  (21)

where G(v) is a piecewise constant function (in units of conductance)and p(v) is also a piecewise constant function (in units of current).The simplest form of an intrinsic conductance model is obtained whenG(v) and p(v) are piecewise constant over just two intervals, as in theHunzinger Cold model described above, for which F(v) takes on thefollowing form:

$\begin{matrix}{{F(v)} = \left\{ \begin{matrix}{\left. {{\frac{C}{\tau_{-}}\left( {v - v_{r}} \right)\mspace{14mu} {if}\mspace{14mu} v} \leq v_{t}}\Rightarrow{G(v)} \right. = {{\frac{C}{\tau_{-}}\mspace{14mu} {and}\mspace{14mu} {p(v)}} = {{- \frac{C}{\tau_{-}}}v_{r}}}} \\{\left. {{\frac{C}{\tau_{+}}\left( {v - v_{t}} \right)\mspace{14mu} {if}\mspace{14mu} v} > v_{t}}\Rightarrow{G(v)} \right. = {{\frac{C}{\tau_{+}}\mspace{14mu} {and}\mspace{14mu} {p(v)}} = {{- \frac{C}{\tau_{+}}}v_{t}}}}\end{matrix} \right.} & (22)\end{matrix}$

All models above, with the exception of the FitzHugh-Nagumo model, arereferred to as two-dimensional hybrid models because, in addition to thedescription given by Eqs. (15) and (16), reset conditions are given.These are provided in such models, for example, because the variable vwill grow to infinity once a voltage threshold is crossed. Therefore thefollowing reset conditions may be used:

If v>v _(peak), then v→v _(r) and u→u+d  (23)

In other words, when the voltage v crosses a spiking threshold v_(peak)(or, for certain aspects, a determination is made that the voltage willcross v_(peak)), the voltage is reset to a resting value v_(r), and therecovery variable u is reset to a value equal to the present value plusa constant d. For certain aspects, u may be reset to a predeterminedconstant value (u_(reset)), instead of u+d. According to certainaspects, a reset condition may occur in response to activation orreception of a control signal. The reset condition in the hybrid modelsenriches the possible behaviors of two-dimensional dynamical systems.

Most generally, the input current I(t) is modeled by a combination ofsynaptic currents I_(SYN)(t) and a generic current I_(ext)(t). Thesynaptic current takes on the form

$\begin{matrix}{{I_{SYN}(t)} = {\sum\limits_{i}{{g_{i}(t)} \cdot \left\lbrack {E_{i} - {v(t)}} \right\rbrack}}} & (24)\end{matrix}$

In Eq. (24), g_(i)(t) indicates the time-dependent conductance for aparticular channel (the i^(th) channel), and E_(i) indicates thereversal potential for that channel. Although a synaptic current of theform above is sufficient to describe substantially linearcurrent-voltage relationships, in some cases (e.g., theN-methyl-D-aspartate (NMDA) channel), the conductance is also a functionof the post-synaptic membrane voltage. In this case Eq. (24) takes on amore complex form, namely

I _(NMDA)(t)=g _(NMDA)(t)·h(v)[E _(NMDA) −v(t)]  (25)

with the function h(v) capturing the dependence on the post-synapticvoltage. The function h(v) is modeled as

h(v)≡1/(1+βe ^(−α·v))  (26)

with the parameters α and β taking on the following values, for example:α=0.062 and β=1/3.57.

Therefore, most generally, up to L different synaptic channels may bemodeled as follows:

$\begin{matrix}{{{I_{SYN}(t)} = {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( {v(t)} \right)}\left( {E_{i} - {v(t)}} \right)}}}{where}} & (27) \\{{h_{i}\left( {v(t)} \right)} = {\frac{1}{1 + {\beta_{i}^{{- \alpha} \cdot {v{(t)}}}}}\mspace{14mu} {with}\mspace{14mu} \left\{ \begin{matrix}{\beta_{i} = \beta} & {{{if}\mspace{14mu} i} = {NMDA}} \\{\beta_{i} = 0} & {otherwise}\end{matrix} \right.}} & (28)\end{matrix}$

Additionally, the time-dependent conductance g_(i)(t) can be modeled bya simple exponential function, an alpha function, or by adifference-of-exponentials function. In the case of a simple decayingexponential with time constant τ where θ(t) is the Heaviside stepfunction, one has

$\begin{matrix}{{g_{i}(t)} = {g_{i}^{- \frac{t}{\tau_{i}}}{\theta (t)}}} & (29)\end{matrix}$

In the case of an alpha function, one has

$\begin{matrix}{{g_{i}(t)} = {g_{i}\frac{t}{\tau_{i}}^{1 - \frac{t}{\tau_{i}}}{\theta (t)}}} & (30)\end{matrix}$

In the case of a difference-of-exponentials function, the exponentialshaving differing rise and decay time constants, one has

$\begin{matrix}{{g_{i}(t)} = {{{\overset{\_}{g}}_{i}\left( {^{- \frac{t}{\tau_{decay}}} - ^{- \frac{t}{\tau_{rise}}}} \right)}{\theta (t)}}} & (31)\end{matrix}$

The constant g _(i) contains a normalization factor such that the peakis equal to g_(i) as follows:

$\begin{matrix}{{{\overset{\_}{g}}_{i} = \frac{g_{i}}{^{- \frac{t_{peak}}{\tau_{decay}}} - ^{- \frac{t_{peak}}{\tau_{rise}}}}}{where}} & (32) \\{t_{peak} = {\frac{\tau_{decay} \cdot \tau_{rise}}{\tau_{decay} - \tau_{rise}} \cdot {\ln \left( \frac{\tau_{decay}}{\tau_{rise}} \right)}}} & (33)\end{matrix}$

In the next sections, approximations to the solution of the dynamicalsystem are explored. Since the system is nonlinear, exact solutionscannot be found. Therefore when referring to or comparing to the “exact”solution (as in Example 1), what is meant is a solution obtainednumerically (e.g., Runge-Kutta), but with great precision.

Formulation and Derivation

Before proceeding with the derivation, some definitions and notationconventions are introduced. Henceforth, boldface capital letters denotematrices, whereas boldface lowercase letters denote vectors. Atwo-dimensional state vector x may be defined as

$\begin{matrix}{x \equiv \begin{bmatrix}v \\u\end{bmatrix}} & (34)\end{matrix}$

The first state variable is the membrane voltage v, and the second statevariable is the recovery variable u. Eqs. (15) and (16) are nowrewritten explicitly for the most general case as follows:

$\begin{matrix}{{C\frac{}{t}{v(t)}} = {{F\left( {v(t)} \right)} - {u(t)} + {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( {v(t)} \right)}\left( {E_{i} - {v(t)}} \right)}} + {I_{ext}(t)}}} & (35) \\{{\frac{}{t}{u(t)}} = {{{abv}(t)} - {a\; {u(t)}} - {abv}_{r}}} & (36)\end{matrix}$

Therefore, in the most general case, one is dealing with a nonlineartime-varying (NLTV) two-dimensional dynamical system, as illustrated at702 in FIG. 7. In the following, such a system is first approximated bya piecewise linear time-varying (pLTV) system at 706 and is successivelyfurther approximated as either: (1) a piecewise linear time-invariantsystem (pLTI) with constant coefficients that are time-intervaldependent at 710 or (2) a pLTI with constant coefficients that aretime-interval independent at 714.

Piecewise Linear Time-Varying System

First the time axis tε[0,∞] is divided into arbitrary non-overlappingtime intervals, namely

{tε[T _(n) ,T _(n+1) ], n=0,1,2,3, . . . }  (37)

Within each time interval (e.g., the n^(th) interval), the neuronbehavior is approximated by a linear time-varying (LTV) system whoseparameters are obtained from the state of the artificial neuron at thestart of the interval (e.g., at time t=T_(n)). In order to develop suchan approximation, first define

$\begin{matrix}{{\Gamma \left( {v,t} \right)} \equiv {\frac{1}{C}\left\{ {{F\left( {v(t)} \right)} + {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( {v(t)} \right)}\left( {E_{i} - {v(t)}} \right)}} + {I_{ext}(t)}} \right\}}} & (38)\end{matrix}$

and the system of Eqs. (35) and (36) becomes

$\begin{matrix}{{\frac{}{t}{v(t)}} = {{\Gamma \left( {v,t} \right)} + {a_{12}{u(t)}}}} & (39) \\{{\frac{}{t}{u(t)}} = {{a_{21}{v(t)}} + {a_{22}{u(t)}} + b_{2}}} & (40)\end{matrix}$

Next, the function Γ(v,t) is approximated with an approximate,affine-linear expression over the interval to tε[T_(n),T_(n+1)]. Thecoefficients in the approximation are a function of the voltage at timeT_(n). To simplify notations when unequivocal, denote v(T_(n))

v_(k) such that

Γ(v,t)≅a ₁₁ [v _(n) ,t]v(t)+b ₁ [v _(n) ,t]  (41)

Some examples of such approximations include the Taylor expansionmethod, the average slope method, the first order linear interpolationmethod, and the optimal linear interpolation method minimizing the L_(p)approximation error. For the Taylor expansion method:

$\begin{matrix}{\left. {{\Gamma \left( {v,t} \right)} \cong {{\Gamma \left( {{v\left( T_{n} \right)},t} \right)} + \frac{\partial{\Gamma \left( {v,t} \right)}}{\partial v}}} \middle| {}_{v{(T_{n})}}\left( {{v(t)} - {v\left( T_{n} \right)}} \right) \right.\left. {{a_{11}\left\lbrack {v_{n},t} \right\rbrack} \equiv \frac{\partial{\Gamma \left( {v,t} \right)}}{\partial v}} \middle| {}_{v_{n}}\mspace{14mu} {and} \right.{{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \equiv {{\Gamma \left( {v_{n},t} \right)} - {{a_{11}\left\lbrack {v_{n},t} \right\rbrack}v_{n}}}}} & \left( {41a} \right)\end{matrix}$

For a given voltage step Δv_(n) in the average slope method, the averageslope is computed as follows:

$\begin{matrix}{{{a_{11}\left\lbrack {{v\left( T_{n} \right)},t} \right\rbrack} \equiv {\frac{{\Gamma \left( {{{v\left( T_{n} \right)} + {\Delta \; v_{n}}},t} \right)} - {\Gamma \left( {{v\left( T_{n} \right)},t} \right)}}{\Delta \; v_{n}}\mspace{14mu} {and}}}{{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \equiv {{\Gamma \left( {v_{n},t} \right)} - {{a_{11}\left\lbrack {v_{n},t} \right\rbrack}v_{n}}}}} & \left( {41b} \right)\end{matrix}$

In the first-order linear interpolation method, the voltage axis ispartitioned into intervals, vε[V_(k),V_(k+1)], k=0, 1, 2, . . . withV_(k)≦v(T_(n))≦V_(k+1). Then

$\begin{matrix}{{{a_{11}\left\lbrack {{v\left( T_{n} \right)},t} \right\rbrack} \equiv {\frac{{\Gamma \left( {V_{k - 1},t} \right)} - {\Gamma \left( {V_{k},t} \right)}}{V_{k + 1} - V_{k}}\mspace{14mu} {and}}}{{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \equiv {{\Gamma \left( {v_{n},t} \right)} - {{a_{11}\left\lbrack {v_{n},t} \right\rbrack}v_{n}}}}} & \left( {41c} \right)\end{matrix}$

In the optimal linear interpolation method minimizing the L_(p)approximation error, the voltage axis is, as before, partitioned intointervals, vε[V_(k), V_(k+1)], k=0, 1, 2, . . . . However, the linearapproximation

{circumflex over (Γ)}(v,t)=a ₁₁ [v _(n) ,t]v(t)+b ₁ [v _(n) ,t]

is designed to minimize a linear approximation error with respect to theoriginal function based on an L_(p) norm. The approximation error overthe interval vε[V_(k),V_(k+1)] may be defined as

${ɛ_{p,n,t}\left( {{a_{11}\left\lbrack {v_{n},t} \right\rbrack},{b_{1}\left\lbrack {v_{n},t} \right\rbrack}} \right)} = \left\lbrack {\int_{V_{k}}^{V_{k + 1}}{{{{\Gamma \left( {v,t} \right)} - {\hat{\Gamma}\left( {v,t} \right)}}}^{p}{v}}} \right\rbrack^{\frac{1}{p}}$where Γ̂(v, t) = a₁₁[v_(n), t]v(t) + b₁[v_(n), t]

Based on Eq. (41), the LTV system described by Eqs. (39) and (40) maynow be expressed as follows:

$\begin{matrix}{{\frac{}{t}{v(t)}} = {{{a_{11}\left\lbrack {v_{n},t} \right\rbrack}{v(t)}} + {a_{12}{u(t)}} + {b_{1}\left\lbrack {v_{n},t} \right\rbrack}}} & (42) \\{{\frac{}{t}{u(t)}} = {{a_{21}{v(t)}} + {a_{22}{u(t)}} + b_{2}}} & (43)\end{matrix}$

More compactly, one has the matrix equation

$\begin{matrix}{{\overset{.}{x}(t)} = {{{A\left( {v_{n},t} \right)}{x(t)}} + {b\left( {v_{n},t} \right)}}} & (44) \\{{\overset{.}{x}(t)} = {{\begin{bmatrix}{a_{11}\left\lbrack {v_{n},t} \right\rbrack} & a_{12} \\a_{21} & a_{22}\end{bmatrix}{x(t)}} + \begin{bmatrix}{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \\b_{2}\end{bmatrix}}} & (45)\end{matrix}$

where a₁₂=−1/C, a₂₁=ab, a₂₂=−a, and b₂=−abv_(r).

The solution of the LTV system above, for tε[T_(n),T_(n+1)], may beexpressed as

$\begin{matrix}{{x(t)} = {{{\Phi \left( {t,T_{n}} \right)}{x\left( T_{n} \right)}} + {\int_{T_{n}}^{t}{{\Phi \left( {\tau,T_{n}} \right)}{b\left( {v_{n},\tau} \right)}{\tau}}}}} & (46)\end{matrix}$

where the transition matrix Φ(t,T_(n)) is given by the Peano-Bakerformula, namely

$\begin{matrix}{{\Phi \left( {t,T_{n}} \right)} = {I + {\int_{T_{n}}^{t}{{A\left( {v_{n},t_{0}} \right)}{t_{0}}}} + {\int_{T_{n}}^{t}{{{A\left( {v_{n},t_{0}} \right)}\left\lbrack {\int_{T_{n}}^{t_{0}}{{A\left( {v_{n},t_{1}} \right)}{t_{1}}}} \right\rbrack}{t_{0}}}} + {\int_{T_{n}}^{t}{{{A\left( {v_{n},t_{0}} \right)}\left\lbrack {\int_{T_{n}}^{t_{0}}{{{A\left( {v_{n},t_{0}} \right)}\left\lbrack {\int_{T_{n}}^{t_{1}}{{A\left( {v_{n},t_{2}} \right)}{t_{2}}}} \right\rbrack}{t_{1}}}} \right\rbrack}{t_{0}}}} + \ldots}} & (47)\end{matrix}$

The series given in Eq. (47) simplifies for some specific cases, forexample:

-   -   (i) A(v_(n),t) is constant and equal to A_(n). Then Φ(t,T_(n))        takes on the familiar form of a matrix exponential, namely

Φ(t,T _(n))=e ^(A) ^(n) ^((t−T) ^(n) ⁾

-   -   -   where the matrix exponential of a matrix A is defined as

$e^{A} = {I + A + \frac{A^{2}}{2!} + \frac{A^{3}}{3!} + \ldots}$

-   -   (ii) A more general case is when A(v_(n),t) and

∫_(T_(n))^(t)A(v_(n), τ)τ

commute for any t.

-   -   -   Then

Φ(t, T_(n)) = exp [∫_(T_(n))^(t)A(v_(n), τ)τ]

-   -   -   Any of the following conditions guarantees that the above            commutative property is satisfied: (a) A(v_(n),t) is            constant; (b) A(v_(n),t)=a(t)M, where a(t) is a scalar            function and M is a constant matrix; or (c)            A(v_(n),t)=Σa_(i)(t)M_(i) where a_(i)(t) is a scalar            function and M_(i)M_(j)=M_(j)M_(i) are constant matrices            that commute for any i,j.

    -   (iii) A somewhat more general case is when for all        tε[T_(n),T_(n+1)], the matrix A(v_(n),t) can be expressed as a        sum of two components:

${A\left( {v_{n},t} \right)} = {{{A_{k}\left( {v_{n},\tau,t_{c}} \right)}\left( {t - t_{c}} \right)^{k}} + {{A_{k + s}\left( {v_{n},\tau,t_{c}} \right)}\left\lbrack \frac{\left( {t - t_{c}} \right)^{k + 1} - \left( {\tau - t_{c}} \right)^{k + 1}}{k + 1} \right\rbrack}^{s}}$

-   -   -   for some integers k≧0 and s≧1, some fixed time-instants            t_(c)≦τ≦t, and some fixed matrices A_(k) (v_(n),τ,t_(c)) and            A_(k) (v_(n),τ,t_(c)). The time-instants and matrices are            “fixed” in the sense that they are independent of the time            parameter tε[T_(n),T_(n+1)]. Perhaps the most useful example            is when k=0 and s=1 or 2. In any case, if the second term in            the above expression (i.e., the term involving the matrix            A_(k+s)(v_(n),τ,t_(c))) is much smaller than the first term            for all tε[T_(n),T_(n+1)], then it turns out that the matrix            Φ(t,T_(n)) may be approximated by the following infinite            series:

${\Phi \left( {t,T_{n}} \right)} = {{\exp \left\lbrack {{A_{k}\left( {v_{n},\tau,t_{c}} \right)}\frac{\; {\left( {t - t_{c}} \right)^{k + 1} - \left( {\tau - t_{c}} \right)^{k + 1}}}{k + 1}} \right\rbrack} + {\sum\limits_{m = 0}^{\infty}{{A_{k}^{m}\left( {v_{n},\tau,t_{c}} \right)}{A_{k + s}\left( {v_{n},\tau,t_{c}} \right)}{G_{m,k,s}\left( {{A_{k}\left( {v_{n},\tau,t_{c}} \right)},\left( {t - t_{c}} \right),\left( {\tau - t_{c}} \right)} \right)}}}}$

-   -   -   where a function G_(m,k,s)(z,t,τ) of a complex variable z            and real-valued parameters t and z is defined. The function            G_(m,k,s)(z,t,τ) is analytic at z=0 and is given by

${G_{m,k,s}\left( {z,t,\tau} \right)} \equiv {\langle{\frac{\partial^{k + s}}{\partial z^{k + s}}\left( {z^{- {({m + 1})}}{\exp \left\lbrack {z\; \frac{t^{k + 1} - \tau^{k + 1}}{k + 1}} \right\rbrack}} \right)}\rangle}$

-   -   -   where the notation <f(z)> denotes the analytical component            of f(z). In the above equation, the analytical component is            obtained by replacing the exponential term with its power            series representation, computing the product and derivative,            and then dropping terms associated with negative powers of            z.

For all other cases of time-varying matrix A(v_(n),t), approximations tothe transition matrix of Eq. (47) can be found, and then approximatesolutions to Eq. (46) may be obtained.

Piecewise Linear Time-Invariant System

Further simplifications may be achieved by transforming the linear timevarying (LTV) system described by Eq. (44) into a linear time-invariant(LTI) system. In order to do so, the coefficient a₁₁[v_(n),t] is heldconstant over the interval tε[T_(n),T_(n+1)]. This result may beachieved in several ways, a few of which are described below. For eachof the affine-linear approximations described above, the constantcoefficient may be defined as follows:

a ₁₁ [v _(n) ,t]≡a ₁₁ [v _(n) ,T _(n)]

for any interval tε[T_(n),T_(n+1)].

Alternatively, if the value of the next time step T_(n+1) is known attime t=T_(n), the average value of the coefficient a₁₁[v_(n),t] may becalculated as follows:

${a_{11}\left\lbrack {v_{n},t} \right\rbrack} \equiv {\frac{1}{T_{n + 1} - T_{n}}{\int_{T_{n}}^{T_{n + 1}}{{a_{11}\left\lbrack {v_{n},\tau} \right\rbrack}{\tau}}}}$

for any interval tε[T_(n),T_(n+1)].

To clarify with an example for the Taylor expansion method, an LTIsystem is obtained by using

${{{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} \equiv \frac{\partial{\Gamma \left( {v,T_{n}} \right)}}{\partial v}}}_{v_{n}}$or$\left. {{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} \equiv {\frac{1}{T_{n + 1} - T_{n}}{\int_{T_{n}}^{T_{n + 1}}\frac{\partial{\Gamma \left( {v,t} \right)}}{\partial v}}}} \middle| {}_{v_{n}}{t} \right.$

The same averaging approach may be applied to any of the methodsdescribed above. The LTI system is now described by the following matrixequation:

$\begin{matrix}{{\overset{.}{x}(t)} = {{\begin{bmatrix}{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} & a_{12} \\a_{21} & a_{22}\end{bmatrix}{x(t)}} + \begin{bmatrix}{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \\b_{2}\end{bmatrix}}} & (48)\end{matrix}$

The solution may be expressed as

$\begin{matrix}{{x(t)} = {{{\Phi \left( {t,T_{n}} \right)}{x\left( T_{n} \right)}} + {\int_{T_{n}}^{t}{{\Phi \left( {{T_{n} + t - \tau_{n}},T_{n}} \right)}{b\left( {v_{n},\tau} \right)}{\tau}}}}} & (49)\end{matrix}$

where the transition matrix Φ(t,T_(n)) is now the matrix exponential

Φ(t,T _(n))=e ^(A(v) ^(n) ^(,T) ^(n) ^()(t−T) ^(n) ⁾  (50)

For the case in which T_(n)=nT (i.e., fixed and uniform time intervalsof length T), Eq. (49) becomes

$\begin{matrix}{{x\left( {{nT} + T} \right)} = {{^{{A{({v_{n},{nT}})}}T}{x({nT})}} + {\int_{nT}^{{nT} + T}{^{{A{({v_{n},{nT}})}}{({{nT} + T - \tau})}}{b\left( {v_{n},\tau} \right)}{\tau}}}}} & (51)\end{matrix}$

In summary, the initial general model is a nonlinear time-varyingtwo-dimensional system. This nonlinear system was first transformed intoa linear time-varying system by applying linear approximations to thenonlinear functional Γ(v, t) over a given time interval. Thetime-varying system was further transformed into a linear time-invariant(LTI) system by taking the time-varying system matrix A(v_(n),t) andapproximating it with a constant matrix A_(n) over the same giveninterval.

Solutions

Now that the problem has been formulated and a number of possible affinelinear approximations to the nonlinear functional have been proposed,the focus now is on the solution to the LTI system over the intervaltε[T_(n),T_(n+1)] given by Eqs. (49) and (51) and their implementation.Further simplifying the notations in Eq. (49) and denotingA(v_(n),T_(n))=A_(n),

$\begin{matrix}{{x(t)} = {{^{A_{n}{({t - T_{n}})}}{x\left( T_{n} \right)}} + {\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}{b\left( {v_{n},\tau} \right)}{\tau}}}}} & (52)\end{matrix}$

and the matrix A(v_(n),T_(n)) and vector b(v_(n),t) are given by

${A_{n} = \begin{bmatrix}{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} & a_{12} \\a_{21} & a_{22}\end{bmatrix}},{{b\left( {v_{n},t} \right)} = \begin{bmatrix}{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \\b_{2}\end{bmatrix}}$

where the coefficients a₁₁[v_(n),T_(n)] and b₁[v_(n),t] are calculatedby any of the methods described above.

Next, the solution for the Taylor expansion method is derived (thesolutions for any of the other methods follow accordingly). In thiscase, one has

$\begin{matrix}{\mspace{20mu} {\left. {{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} \equiv \frac{\partial{\Gamma \left( {v,T_{n}} \right)}}{\partial v}} \middle| {}_{v_{n}}\mspace{14mu} {and} \right.\mspace{20mu} {{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \equiv {{\Gamma \left( {v_{n},t} \right)} - {{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack}v_{n}}}}\mspace{20mu} {and}}} & (53) \\{\mspace{20mu} {{\Gamma \left( {v_{n},t} \right)} \equiv {\frac{1}{C}\left\{ {{F\left( v_{n} \right)} + {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)}} + {I_{{ex}\; t}(t)}} \right\}}}} & (54) \\{\left. {{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} \equiv \frac{\partial{\Gamma \left( {v,T_{n}} \right)}}{\partial v}} \right|_{v_{n}} = \left. {\frac{1}{C}\frac{\partial{F(v)}}{\partial v}} \middle| {}_{v_{n}}{{+ \frac{1}{C}}{\sum\limits_{i = 1}^{L}{{g_{i}\left( T_{n} \right)}\left\lbrack \frac{\partial{h_{i}(v)}}{\partial v} \middle| {}_{v_{n}}{\left( {E_{i} - v_{n}} \right) - {h_{i}\left( v_{n} \right)}} \right\rbrack}}} \right.} & (55)\end{matrix}$

Simplifying notations,

$\begin{matrix}{{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack} = {{\frac{1}{C}{F^{\prime}\left( v_{n} \right)}} + {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{g_{i}\left( T_{n} \right)}\left\lbrack {{{h_{i}^{\prime}\left( v_{n} \right)}\left( {E_{i} - V_{n}} \right)} - {h_{i}\left( v_{n} \right)}} \right\rbrack}}}}} & (56)\end{matrix}$

It is convenient to express the vector b(v_(n),t) as follows:

${b\left( {v_{n},t} \right)} = {\begin{bmatrix}{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \\b_{2}\end{bmatrix} = {\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix} + \begin{bmatrix}{b_{1}(t)} \\0\end{bmatrix} + \begin{bmatrix}{\frac{1}{C}{I_{ext}(t)}} \\0\end{bmatrix}}}$ ${{where}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}} = {{\begin{bmatrix}{\frac{1}{C}\left\{ {{F\left( v_{n} \right)} - {{a_{11}\left\lbrack {v_{n},T_{n}} \right\rbrack}v_{n}}} \right\}} \\b_{2}\end{bmatrix}\mspace{14mu} {{and}\text{}\begin{bmatrix}{b_{1}(t)} \\0\end{bmatrix}}} = \begin{bmatrix}{\frac{1}{C}\left\{ {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)}} \right\}} \\0\end{bmatrix}}$

Eq. (52) may then be expressed as

$\begin{matrix}{{x(t)} = {{^{A_{n}{({t - T_{n}})}}{x\left( T_{n} \right)}} + {\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}{{\tau \begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}}}}} + {\int_{T_{n}}^{t}{{^{A_{n}{({t - \tau})}}\begin{bmatrix}{b_{1}(\tau)} \\0\end{bmatrix}}{\tau}}} + {\int_{T_{n}}^{t}{{^{A_{n}{({t - \tau})}}\begin{bmatrix}{\frac{1}{C}I_{ext}(\tau)} \\0\end{bmatrix}}{\tau}}}}} & (57)\end{matrix}$

The first integral of Eq. (57) may be solved in closed form, such thatone obtains

$\begin{matrix}{{\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}{{\tau \begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}}}}} = {{A_{n}^{- 1}\left( {^{A_{n}{({t - T_{n}})}} - I} \right)}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}}} & (58)\end{matrix}$

where I is the 2×2 identity matrix.

Eq. (52) may now be expressed as follows:

$\begin{matrix}{{x(t)} = {{^{A_{n}{({t - T_{n}})}}{x\left( T_{n} \right)}} + {{A_{n}^{- 1}\left( {^{A_{n}{({t - T_{n}})}} - I} \right)}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}} + {\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}{b_{1}(\tau)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}} + {\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}\frac{1}{C}{I_{ext}(\tau)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}} & (59)\end{matrix}$

If the external current is such that closed-form expression for the lastintegral of Eq. (59) cannot be obtained, it may have to be approximated.One example of an approximation is the zero-hold, namely

$\begin{matrix}{{{\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}\frac{1}{C}{I_{ext}(\tau)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}} \cong {\frac{1}{C}{I_{ext}\left( T_{n} \right)}{\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}} = {\frac{1}{C}{I_{ext}\left( T_{n} \right)}{{A_{n}^{- 1}\left( {^{A_{n}{({t - T_{n}})}} - I} \right)}\begin{bmatrix}1 \\0\end{bmatrix}}}} & (60)\end{matrix}$

In this case, the closed-form solution for any intervaltε[T_(n),T_(n+1)] takes on the following form:

$\begin{matrix}{{x(t)} = {{^{A_{n}{({t - T_{n}})}}{x\left( T_{n} \right)}} + {{A_{n}^{- 1}\left( {^{A_{n}{({t - T_{n}})}} - I} \right)}\begin{bmatrix}{b_{0} + {\frac{1}{C}{I_{ext}\left( T_{n} \right)}}} \\b_{2}\end{bmatrix}} + {\int_{T_{n}}^{t}{^{A_{n}{({t - \tau})}}{b_{1}(\tau)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}} & (61)\end{matrix}$

A key ingredient in Eq. (61) is the matrix exponential e^(A) ^(n) ^(t)which can be expressed as follows. Let λ_(n1) and λ_(n2) be theeigenvalues of the 2×2 matrix A_(n). Then

$\begin{matrix}{\mspace{79mu} {^{A_{n}t} = {{{^{\lambda_{n\; 1}t}I} + {\frac{^{\lambda_{n\; 1}t} - ^{\lambda_{n\; 2}t}}{\lambda_{n\; 1} - \lambda_{n\; 2}}\left( {A_{n} - {\lambda_{n\; 1}I}} \right)\mspace{14mu} {for}\mspace{14mu} \lambda_{n\; 1}}} \neq {\lambda_{n\; 2}\mspace{14mu} {real}}}}} & (62) \\{\mspace{79mu} {^{A_{n}t} = {{{^{\lambda_{n\; 1}t}I} + {t\; {^{\lambda_{n\; 1}t}\left( {A_{n} - {\lambda_{n\; 1}I}} \right)}\mspace{14mu} {for}\mspace{14mu} \lambda_{n\; 1}}} = {\lambda_{n\; 2}\mspace{14mu} {real}}}}} & (63) \\{^{A_{n}t} = {{{^{\alpha \; t}\cos \; \beta \; t\; I} + {\frac{^{\alpha \; t}\sin \; \beta \; t}{\beta}\left( {A_{n} - {\alpha \; I}} \right)\mspace{14mu} {for}\mspace{14mu} \lambda_{n\; 1}}} = {\lambda_{n\; 2}^{*} = {{\alpha + {\; \beta \mspace{14mu} \beta}} > 0}}}} & (64)\end{matrix}$

where λ_(n1) and λ_(n2) are the solutions of the characteristic equationdet(A_(n)−λI)=0.

Eqs. (62)-(64) give explicit expressions for the matrix exponential.Alternatively, a large number of efficient numerical algorithms areavailable to compute the matrix exponential. Therefore, it is feasibleto determine the state of the system at time t=T_(n+1) through eithermethod. But in a hardware implementation, it may be beneficial to use adifferent approach, as described next.

First, notice that in the matrix A_(n), only one coefficient a₁₁[v_(n),T_(n)] changes as the system evolves; the remaining threecoefficients remain constant throughout. The coefficient is given, forexample, by Eq. (56). An example efficient procedure is as follows: (1)quantize the first state variable (i.e., the membrane voltage v_(n)) ateach time to step to a finite number of values and (2) index lookuptables with pre-calculated values of the functions F′(v_(n)h_(i)′(v_(n)), and h_(i)(v_(n)). Similarly, time instants may bequantized, as well, and pre-calculated values of g_(i)(T_(n)) may bestored in lookup tables. The retrieved values may now be used toevaluate a₁₁[v_(n),T_(n)] via Eq. (56) and to quantize the result to afinite number of values Q(a₁₁[v_(n),T_(n)]).

Similarly, a number of choices for ΔT_(n) may be pre-selected, and thematrix exponential and the matrix inverse of the matrices appearing inEq. (62) may be pre-calculated and addressed by Q(a₁₁[v_(n),T_(n)]) andΔT_(n).

Adaptive and Fixed Time Steps

In the previous section, the expressions for the state of the LTI systemwere derived at arbitrary time instants. The following describes in moredetail the case where the time steps are chosen adaptively and thesimpler case of a synchronous system where the time instants are uniformwith a fixed and predetermined step size ΔT_(n)=T.

An adaptive strategy may be beneficial both with respect to complexityand precision. For example, time steps may be chosen farther apart(i.e., larger ΔT_(n)) when the state is evolving slowly, thus increasingthe speed of computations. In contrast, the time steps may be made smallwhen, for example, the neuron model is near a spiking event, therebyincreasing the precision with which the spike time is calculated. Next,a possible algorithm for the choice of the time steps is described.

As Eq. (63) shows, the time constants of the system evolution areinversely proportional to the eigenvalues of the matrices A_(n). Inparticular, if the eigenvalues are close to zero, then the timeconstants are large, and the system evolves slowly. If, on the contrary,at least one of the eigenvalues becomes large (and possibly positive),the system evolves very rapidly. One example algorithm for the choice ofthe time steps that accomplishes the above task is the following:

$\begin{matrix}{T_{n + 1} = {T + {\min \left( {\frac{\ln \; \mu_{n}}{\lambda_{n,\max}},{\Delta \; T_{\max}}} \right)}}} & (65)\end{matrix}$

In Eq. (65), λ_(n,max) denotes the maximum eigenvalue (in magnitude) ofthe matrix A_(n), while ΔT_(max) and μ_(n) are configurable parameters.This ensures that the following relationships hold:

λ_(n,max) ≦|a ₁₁ [v(T _(n)),T _(n) ]|+|a ₂₂ |, T _(n+1) −T _(n) ≦ΔT_(max)  (66)

and for any interval tε[T_(n), T_(n+1)]

∥x(t)∥≦μ_(n) ∥x(T _(n))∥  (67)

Next, the general solutions of the previous section are specialized forthe important case of fixed and uniform time steps of size T. Startingwith Eq. (61), one obtains

$\begin{matrix}{{x\left( {{n\; T} + T} \right)} = {{^{A_{n}T}{x\left( {n\; T} \right)}} + {{A_{n}^{- 1}\left( {^{A_{n}T} - I} \right)}\begin{bmatrix}{b_{0} + {\frac{1}{C}{I_{ext}\left( {n\; T} \right)}}} \\b_{2}\end{bmatrix}} + {\int_{0}^{T}{^{A_{n}t}{b_{1}\left( {{n\; T} + T - t} \right)}{{t\begin{bmatrix}1 \\0\end{bmatrix}}}}}}} & (68)\end{matrix}$

Example 1 Taylor Expansion for the Izhikevich (Simple) Model

In this example, Eq. (68) is evaluated for the simple model andlinearization based on the Taylor expansion method. To simplify theformulas somewhat, it is assumed that the external current is notpresent.

First, the coefficients of the Taylor expansion for the simple model,i.e., F(v)=k(v−v_(r))(v−v_(t)), are derived starting from Eq. (56).Provided that

$\begin{matrix}{\left. {\frac{1}{C}\frac{\partial{F(v)}}{\partial v}} \right|_{v_{n}} = {{\frac{2\; k}{C}v_{n}} - {\frac{k}{C}\left( {v_{r} + v_{t}} \right)}}} & (69)\end{matrix}$

one obtains

$\begin{matrix}{{a_{11}\left\lbrack {v_{n},{n\; T}} \right\rbrack} = {{\frac{2\; k}{C}v_{n}} - {\frac{k}{C}\left( {v_{r} + v_{t}} \right)} + {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{g_{i}\left( {n\; T} \right)}\left\lbrack {{{h_{i}^{\prime}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)} - {h_{i}\left( v_{n} \right)}} \right\rbrack}}}}} & (70) \\{{b_{1}\left\lbrack {v_{n},t} \right\rbrack} = {{{\Gamma \left( {v_{n},t} \right)} - {{a_{11}\left\lbrack {v_{n},{n\; T}} \right\rbrack}v_{n}}} = \begin{matrix}{{\frac{1}{C}\left\{ {{{k\left( {v_{n} - v_{r}} \right)}\left( {v_{n} - v_{t}} \right)} + {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)}}} \right\}} -} \\\left\{ {{\frac{2\; k}{C}v_{n}^{2}} - {\frac{k}{C}\left( {v_{r} + v_{t}} \right)v_{n}} + {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{{g_{i}\left( {n\; T} \right)}\left\lbrack {{{h_{i}^{\prime}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)} - {h_{i}\left( v_{n} \right)}} \right\rbrack}v_{n}}}}} \right\}\end{matrix}}} & (71)\end{matrix}$

For non-NMDA synaptic channels, Eqs. (70) and (71) simplify to

$\begin{matrix}{\mspace{79mu} {{a_{11}\left\lbrack {v_{n},{n\; T}} \right\rbrack} = {{\frac{2k}{C}v_{n}} - {\frac{k}{C}\left( {v_{r} + v_{t}} \right)} - {\frac{1}{C}{\sum\limits_{i = 1}^{L}{g_{i}\left( {n\; T} \right)}}}}}} & (72) \\{{b_{1}\left\lbrack {v_{n},{n\; T}} \right\rbrack} = {{\frac{1}{C}\left\{ {{{- k}\; v_{n}^{2}} + {k\; v_{r}v_{t}} + {\sum\limits_{i = 1}^{L}{{g_{i}\left( {n\; T} \right)}v_{n}}}} \right\}} + {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{g_{i}(t)}\left( {E_{i} - v_{n}} \right)}}}}} & (73)\end{matrix}$

Taking

$b_{0} = {\frac{1}{C}\left\{ {{{- k}\; v_{n}^{2}} + {k\; v_{r}v_{t}} + {\sum\limits_{i = 1}^{L}{{{g_{i}\left( {n\; T} \right)}\left\lbrack {{{h_{i}^{\prime}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)} - {h_{i}\left( v_{n} \right)}} \right\rbrack}v_{n}}}} \right\}}$

Eq. (68) now becomes

$\begin{matrix}{{{x\left( {{n\; T} + T} \right)} = {{^{A_{n}T}{x\left( {n\; T} \right)}} + {{A_{n}^{- 1}\left( {^{A_{n}T} - I} \right)}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}} + {\sum\limits_{i = 1}^{L}{\left\lbrack {{h_{i}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)} \right\rbrack {\int_{0}^{T}{^{A_{n}\tau}{g_{i}\left( {{n\; T} + T - \tau} \right)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}}}}\mspace{20mu} {with}{A_{n} = {\quad{{\begin{bmatrix}{{\frac{2k}{C}v_{n}} - {\frac{k}{C}\left( {v_{r} + v_{t}} \right)} + {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{g_{i}\left( {n\; T} \right)}\begin{bmatrix}{{{h_{i}^{\prime}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)} -} \\{h_{i}\left( v_{n} \right)}\end{bmatrix}}}}} & {- \frac{1}{C}} \\{ab} & {- a}\end{bmatrix}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}} = {\quad\begin{bmatrix}{\frac{1}{C}\left\{ {{{- k}\; v_{n}^{2}} + {k\; v_{r}v_{t}} + {\sum\limits_{i = 1}^{L}{{{g_{i}\left( {n\; T} \right)}\left\lbrack {{{h_{i}^{\prime}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)} - {h_{i}\left( v_{n} \right)}} \right\rbrack}v_{n}}}} \right\}} \\{- {abv}_{r}}\end{bmatrix}}}}}} & (74)\end{matrix}$

Once again for non-NMDA synaptic currents, Eq. (74) simplifies to

$\begin{matrix}{{{x\left( {{n\; T} + T} \right)} = {{^{A_{n}T}{x\left( {n\; T} \right)}} + {{A_{n}^{- 1}\left( {^{A_{n}T} - I} \right)}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}} + {\sum\limits_{i = 1}^{L}{E_{i}{\int_{0}^{T}{^{A_{n}\tau}{g_{i}\left( {{n\; T} + T - \tau} \right)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}}}}\mspace{20mu} {with}\mspace{20mu} {A_{n} = {{\begin{bmatrix}{{\frac{2k}{C}v_{n}} - {\frac{k}{C}\left( {v_{r} + v_{t}} \right)} - {\frac{1}{C}{\sum\limits_{i = 1}^{L}{g_{i}\left( {n\; T} \right)}}}} & {- \frac{1}{C}} \\{ab} & {- a}\end{bmatrix}\mspace{20mu}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}} = {\quad\begin{bmatrix}{\frac{1}{C}\left\{ {{{- k}\; v_{n}^{2}} + {k\; v_{r}v_{t}} + {\sum\limits_{i = 1}^{L}{{g_{i}\left( {n\; T} \right)}v_{n}}}} \right\}} \\{- {abv}_{r}}\end{bmatrix}}}}} & (75)\end{matrix}$

Furthermore, the integrals in Eqs. (74) and (75) can be solved in closedform. For example, if

$\begin{matrix}{\mspace{79mu} {{{g_{i}(t)} = {g_{i}^{{- t}/\tau_{i}}{\theta (t)}}}\mspace{20mu} {then}{{\int_{0}^{T}{^{A_{n}\tau}{g_{i}\left( {{n\; T}\; + T - \tau} \right)}{\tau}}} = {{g_{i}\left( {{n\; T} + T} \right)}\left( {^{A_{n}} + {\tau_{i}^{- 1}I}} \right)^{- 1}\left( {^{{({A_{n} + {\tau_{i}^{- 1}I}})}T} - I} \right)}}}} & (76)\end{matrix}$

To complete the example, the accuracy of Eq. (75) is tested with T=1 msfor a slow excitatory neuron model with parameters k=0.7, C=100, a=0.03,b=−2, v_(r)=−60 mV, and v_(t)=−40 mV and a single exponential AMPAchannel with parameters τ_(AMPA)=5 ms, E_(AMPA)=0 mV, and g_(AMPA)=5. Asillustrated in the membrane potential plot 500 and the recovery currentplot 520 of FIGS. 5A and 5B, respectively, the linearization based onEq. (75) accurately tracks the plots based on the nonlinear time-varyingmodel.

Example 2 Subthreshold Dynamics of the Hunzinger Cold Model

Another example is developed in an effort to examine the subthresholddynamics of the simple intrinsic conductance model known as theHunzinger Cold model. In this example, no synaptic currents and thesimple, but interesting case of impulsive external current are assumed.

In the Hunzinger Cold model (see Eq. (22) above) when the membranevoltage is below threshold, the matrix, A_(n) is constant and equal toA⁻:

$A_{-} \equiv \begin{bmatrix}{{- 1}/\tau_{-}} & {{- 1}/C} \\{ab} & {- a}\end{bmatrix}$

The derivation can be further simplified if the first state variable isdefined as the membrane voltage minus the reference voltage. With such adefinition, all the constant terms are equal to zero, and Eq. (61)simplifies to

$\begin{matrix}{{x(t)} = {{^{A_{-}{({t - T_{n}})}}{x\left( T_{n} \right)}} + {\frac{1}{C}{\int_{T_{n}}^{t}{^{A_{-}{({t - \tau})}}{I_{ext}(\tau)}{{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}}} & (77)\end{matrix}$

Furthermore, if the external current is assumed to be a Dirac deltafunction at time T_(n) with amplitude I, i.e., I_(ext)=(t)+Iδ(t−T_(n)),then

$\begin{matrix}{{x(t)} = {^{A_{-}{({t - T_{n}})}}\left( {{x\left( T_{n} \right)} + {\frac{1}{C}\begin{bmatrix}1 \\0\end{bmatrix}}} \right)}} & (78)\end{matrix}$

Note that impulsive inputs have the same effect on the system as that ofan initial condition. If the eigenvalues are real and λ₁≠λ₂, Eq. (62)may be substituted in Eq. (78) to obtain the following continuous timesolution for v<v_(t):

$\begin{matrix}{{x(t)} = {\left\lbrack {{^{\lambda_{1}{({t - T_{n}})}}I} + {\frac{^{\lambda_{1}{({t - T_{n}})}} - ^{\lambda_{2}{({t - T_{n}})}}}{\lambda_{1} - \lambda_{2}}\left( {A_{-} - {\lambda_{1}I}} \right)}} \right\rbrack \left( {{x\left( T_{n} \right)} + {\frac{I}{C}\begin{bmatrix}1 \\0\end{bmatrix}}} \right)}} & (79)\end{matrix}$

By defining a new matrix

${\overset{\_}{A}}_{-} = \frac{\left( {A_{-} - {\lambda_{1}I}} \right)}{\lambda_{1} - \lambda_{2}}$

and denoting

${{x\left( T_{n} \right)} + {\frac{I}{C}\begin{bmatrix}1 \\0\end{bmatrix}}} = \begin{bmatrix}v_{n} \\u_{n}\end{bmatrix}$

Eq. (72) may be written as

$\begin{matrix}{{x(t)} = {\left\lbrack {{^{\lambda_{1}{({t - T_{n}})}}\left( {I + {\overset{\_}{A}}_{-}} \right)} - {^{\lambda_{2}{({t - T_{n}})}}{\overset{\_}{A}}_{-}}} \right\rbrack \begin{bmatrix}v_{n} \\u_{n}\end{bmatrix}}} & (80)\end{matrix}$

In particular for T_(n)=0, one obtains

$\begin{matrix}{{x(t)} = {\left\lbrack {{^{\lambda_{1}t}\left( {I + {\overset{\_}{A}}_{-}} \right)} - {^{\lambda_{2}t}{\overset{\_}{A}}_{-}}} \right\rbrack \begin{bmatrix}{v_{0} + {I/C}} \\u_{0}\end{bmatrix}}} & (81)\end{matrix}$

The membrane potential plot 600 and the recovery current plot 620 ofFIGS. 6A and 6B, respectively, show the time evolution of Eq. (81) forthe following values of the parameters: τ⁻=14 ms, C=100, a=0.03, b=−2,I_(ext)(t)=100δ(t), and v₀=u₀=0. The membrane voltage and recoverycurrent are given by the following functions:

v(t)−v _(r)=0.8228e ^(−0.0827t)+0.1771e ^(−0.0186t)  (82)

u(t)=0.9351e ^(−0.08277t)+0.9351e ^(−0.0186t)  (83)

Since the matrices and eigenvalues involved in the computation of Eq.(80) are known a priori, the time course of the system may be calculatedvery easily and with the desired precision by pre-calculating andstoring (e.g., in a look-up table) the two exponential functions with agiven time resolution. Alternatively, fast algorithms for thecomputation (approximation) of the exponential function may be employed,and Eq. (80) may be computed in real time.

The above formulation is particularly relevant in an event-basedimplementation where it is of interest to compute the state vector atthe time of an event. For instance given Eq. (80) and the occurrence ofan event at time T_(n+1), one may wish to compute T_(n+1) ⁻. WithΔ_(n)≡T_(n+1)−T_(n), one has

$\begin{matrix}{{x\left( T_{n + 1}^{-} \right)} = {\left\lbrack {{^{\lambda_{1}\Delta_{n}}\left( {I + {\overset{\_}{A}}_{-}} \right)} - {^{\lambda_{2}\Delta_{n}}{\overset{\_}{A}}_{-}}} \right\rbrack \begin{bmatrix}v_{n} \\u_{n}\end{bmatrix}}} & (84)\end{matrix}$

One important event is the spike time (i.e. the time at which theartificial neuron's membrane voltage crosses a peak threshold, asdescribed above). When the model for the artificial neuron is close tothe spike time, at least one of the eigenvalues is positive anddominant. Therefore, a good approximation may be achieved by retainingonly the dominant term in the expression of the voltage. In other words,since the function v(t) is of the form shown in Eq. (82), i.e.,

v(t)=c ₁ e ^(λ) ¹ ^(t) +c ₂ e ^(λ) ² ^(t)  (85)

v(t) may be approximated by its dominant term, i.e.,

v(t)≅c ₁ e ^(λ) ¹ ^(t) +c ₂  (86)

The spiking time may be easily obtained from

$\begin{matrix}{t_{spike} \cong {\frac{1}{\lambda_{1}}{\ln \left( {\left( {v_{peak} - c_{2}} \right)/c_{1}} \right)}}} & (87)\end{matrix}$

Additional Simplifications

The end of the subsection entitled “Solutions” mentions that in order topre-calculate the matrix exponential and matrix inverse, the firstelement a₁₁[v_(n),T_(n)] of the matrix A_(n) should be quantized to afinite number of values that cover the dynamic range. As seen from thederivations of the previous subsection and specifically from Example 1,the coefficient a₁₁[v_(n),T_(n)] is a function of v_(n) and of theconductance at time T_(n). Although the dynamic range of the membranevoltage is fairly limited, the dynamic range of the conductance can bevery large. Therefore, a large number of pre-calculated matrices maymost likely have to be stored in this case. Next, furthersimplifications are explored in order to reduce the memory demands. Oneapproach is that of pre-computing and storing the desired matrices overa limited range and performing calculations to approximate the desiredmatrices for values beyond that range.

For example, when the coefficient a₁[v_(n),T_(n)] becomes much larger(in magnitude) than the other three fixed coefficients of the matrixA_(n), the following approximation may be used:

$\begin{matrix}{{\psi (A)} \cong {{{\frac{1}{\left( {a_{11} - a_{22}} \right)^{2} + {a_{12}a_{21}}}\begin{bmatrix}{a_{11} - a_{22}} & {- a_{12}} \\a_{21} & {a_{11} - a_{22}}\end{bmatrix}}\begin{bmatrix}{\psi \left( \gamma_{1} \right)} & 0 \\0 & {\psi \left( \gamma_{2} \right)}\end{bmatrix}}{\quad{{\begin{bmatrix}{a_{11} - a_{22}} & a_{12} \\{- a_{21}} & {a_{11} - a_{22}}\end{bmatrix}\mspace{20mu} \gamma_{1}} = {{a_{11} + {\frac{a_{12}a_{21}}{a_{11} - a_{22}}\mspace{14mu} {and}\mspace{14mu} \gamma_{2}}} = {a_{22} - \frac{a_{12}a_{21}}{a_{11} - a_{22}}}}}}}} & (88)\end{matrix}$

for any complex-valued function ψ(·), such as ψ(x)=e^(x) and ψ(x)=x⁻¹.With this approach, memory demands and computational complexity can betraded off.

Further reduction of the memory demands may be accomplished via certainapproximations. For example, an affine linear approximation may beapplied to a revised functional

${\Gamma \left( {v,t} \right)} = {{\Gamma (v)} \equiv {\frac{1}{C}{F\left( {v(t)} \right)}}}$

rather than that defined by Eq. (38). Additionally, the synaptic currentmay be treated as an external current, and the resulting integral may beapproximated by a backward rectangular rule. In this case, the dynamicrange involved for pre-calculating and storing the matrix exponentialand matrix inverse is that of the membrane voltage alone and is,therefore, quite small. The tradeoff is one of memory demands versusaccuracy of the solutions.

Next, the paragraph above is explained in more detail by re-deriving thesolution under this additional simplification. One may begin as beforewith Eqs. (35) and (36), repeated here for convenience:

$\begin{matrix}{{C\frac{}{t}{v(t)}} = {{F\left( {v(t)} \right)} - {u(t)} + {\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( {v(t)} \right)}\left( {E_{i} - {v(t)}} \right)}} + {I_{{ext}\;}(t)}}} & (89) \\{{\frac{}{t}{u(t)}} = {{{abv}(t)} - {a\; {u(t)}} - {abv}_{r}}} & (90)\end{matrix}$

Now, however, the functional in Eq. (38) is modified as follows:

$\begin{matrix}{{\Gamma \left( {v,t} \right)} = {{\Gamma (v)} \equiv {\frac{1}{C}{F\left( {v(t)} \right)}}}} & (91)\end{matrix}$

and the system of equations is reformulated as

$\begin{matrix}{{\frac{\;}{t}{v(t)}} = {{\Gamma (v)} + {a_{12}{u(t)}} + {\frac{1}{C}\left\{ {{\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( {v(t)} \right)}\left( {E_{i} - {v(t)}} \right)}} + {I_{ext}(t)}} \right\}}}} & (92) \\{\mspace{85mu} {{\frac{\;}{t}{u(t)}} = {{a_{21}{v(t)}} + {a_{22}{u(t)}} + b_{2}}}} & (93)\end{matrix}$

The system matrix A_(n) now simplifies to

$\begin{matrix}{A_{n} \equiv \begin{bmatrix}{a_{11}\left\lbrack v_{n} \right\rbrack} & a_{12} \\a_{21} & a_{22}\end{bmatrix}} & (94)\end{matrix}$

with the vector b(v_(n),t) is now defined as

$\begin{matrix}\begin{matrix}{{b\left( {v_{n},t} \right)} = \begin{bmatrix}{b_{1}\left\lbrack {v_{n},t} \right\rbrack} \\b_{2}\end{bmatrix}} \\{= \begin{bmatrix}{\frac{1}{C}\left\{ {{\sum\limits_{i = 1}^{L}{{g_{i}(t)}{h_{i}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)}} + {I_{ext}(t)}} \right\}} \\b_{2}\end{bmatrix}}\end{matrix} & (95)\end{matrix}$

With the above approximation, the system matrix (coefficient a₁₁)dependence on the conductance values is avoided, thus considerablyreducing the dynamic range to be covered in computing the matrixexponential. The synaptic current may be treated as an external currentby holding the membrane voltage to a constant value for the time period.In a fixed-time-step implementation, the solution would then be

x(nT+T)=e ^(A) ^(n) ^(T) x(nT)+Σ_(nT) ^(nT+T) e ^(A) ^(n) ^((nT+T−τ))b(v _(n),τ)dτ  (96)

The first component of the vector integral in Eq. (96) is

$\begin{matrix}{\int_{nT}^{{nT} - T}{^{A_{n}{({{nT} + T - \tau})}}\frac{1}{C}\left\{ {{\sum\limits_{i = 1}^{L}{{g_{i}(\tau)}{h_{i}\left( v_{n} \right)}\left( {E_{i} - v_{n}} \right)}} + {I_{ext}(t)}} \right\} \ {\tau}}} & (97)\end{matrix}$

which can be solved in most cases of interest. The best results for theabove approximation have been obtained by using a backward rectangularrule, namely using the value v_(n)=v(nT+T) in Eq. (97).

Example 3 Approximation with a Backward Rectangular Rule

To simplify the exposition and notations, this example assumes that thesynaptic current does not contain voltage-dependent conductance channels(e.g., NMDA channels) and that there is no external current I_(ext)(t).The derivation is obtained for fixed step sizes of length T. To begin,Eq. (68) may be simplified as follows:

$\begin{matrix}{{x\left( {{nT} + T} \right)} = {{^{A_{n}T}{x({nT})}} + q + {\frac{1}{C}{\int\limits_{0}^{T}{^{A_{n^{\tau}}}{\sum\limits_{i = 1}^{L}{{g_{i}\left( {{nT} + T - \tau} \right)}\left( {E_{i} - {v\left( {{nT} + T} \right)}} \right){{\tau \begin{bmatrix}1 \\0\end{bmatrix}}}}}}}}}} & (98)\end{matrix}$

where the vector q contains constant terms, namely

$\begin{matrix}{q = {{A_{n}^{- 1}\left( {^{A_{n}T} - I} \right)}\begin{bmatrix}b_{0} \\b_{2}\end{bmatrix}}} & (99)\end{matrix}$

For notational convenience, the following vector and matrix may also bedefined:

$\begin{matrix}{{b \equiv \begin{bmatrix}1 \\0\end{bmatrix}}{B \equiv \begin{bmatrix}1 & 0 \\0 & 0\end{bmatrix}}} & (100)\end{matrix}$

such that Eq. (98) may be rewritten as

$\begin{matrix}{{x\left( {{nT} + T} \right)} = {{^{A_{n}T}{x({nT})}} + {\frac{1}{C}{\int\limits_{0}^{T}{^{A_{n}T}{\sum\limits_{i = 1}^{L}{{g_{i}\left( {{nT} + T - \tau} \right)}E_{i}{\tau}\; b}}}}} - {\frac{1}{C}{\int\limits_{0}^{T}{^{A_{n}\tau}{\sum\limits_{i = 1}^{L}{{g_{i}\left( {{nT} + T - \tau} \right)}{\tau}\; {B \cdot {x\left( {{nT} + T} \right)}}}}}}} + q}} & (101)\end{matrix}$

Now, the solution is derived for g_(i)(t) modeled as simple exponentialsper Eq. (29). In this case, one has

$\begin{matrix}{{x\left( {{nT} + T} \right)} = {{^{A_{n}T}{x({nT})}} + {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{g_{i}({nT})}E_{i}{\int\limits_{0}^{T}{^{A_{n}\tau}^{{- {({T - \tau})}}/\tau_{i}}{{\tau} \cdot b}}}}}} - {\frac{1}{C}{\sum\limits_{i = 1}^{L}{{g_{i}({nT})}{\int\limits_{0}^{T}{^{A_{n}\tau}^{{- {({T - \tau})}}/\tau_{i}}{\tau}\; {B \cdot {x\left( {{nT} + T} \right)}}}}}}} + q}} & (102)\end{matrix}$

The integrals can be readily solved, and if one further defines

$\begin{matrix}{{H_{i} \equiv {\int\limits_{0}^{T}{^{A \cdot \sigma}^{\frac{\sigma - T}{\tau_{i}}}}}}{{{d\; \sigma} = \begin{bmatrix}h_{11}^{i} & h_{12}^{i} \\h_{21}^{i} & h_{22}^{i}\end{bmatrix}},{{h_{i} \equiv {H_{i} \cdot b}} = \begin{bmatrix}h_{11}^{i} \\h_{21}^{i}\end{bmatrix}}}} & (103)\end{matrix}$

then Eq. (102) may be written as

$\begin{matrix}{{x\left( {{nT} + T} \right)} = {{^{A \cdot T}{x({nT})}} + {\sum\limits_{i}^{\;}{\frac{E_{i}{g_{i}({nT})}}{C} \cdot h_{i}}} - {\sum\limits_{i}^{\;}{\frac{g_{i}({nT})}{C} \cdot H_{i} \cdot B \cdot {x\left( {{nT} + T} \right)}}} + q}} & (104)\end{matrix}$

Collecting terms containing the state at time nT+T on the left-hand sideleads to

$\begin{matrix}{{\left\lbrack {I + {\sum\limits_{i}^{\;}{\frac{g_{i}({nT})}{C} \cdot H_{i} \cdot B}}} \right\rbrack \cdot {x\left( {{nT} + T} \right)}} = {{^{A \cdot T}{x({nT})}{\sum\limits_{i}^{\;}{\frac{E_{i}{g_{i}({nT})}}{C} \cdot h_{i}}}} + q}} & (105)\end{matrix}$

and finally

$\begin{matrix}{{\cdot {x\left( {{nT} + T} \right)}} = {\left\lbrack {I + {\sum\limits_{i}^{\;}{\frac{g_{i}({nT})}{C} \cdot H_{i} \cdot B}}} \right\rbrack^{- 1} \cdot \left\{ {{^{A \cdot T}{x({nT})}} + {\sum\limits_{i}^{\;}{\frac{E_{i}{g_{i}({nT})}}{C} \cdot h_{i}}} + q} \right\}}} & (106)\end{matrix}$

Notice that the matrix to be inverted is

$\begin{matrix}\begin{matrix}{\left\lbrack {I + {\sum\limits_{i}^{\;}{\frac{g_{i}({nT})}{C} \cdot H_{i} \cdot B}}} \right\rbrack = {\begin{bmatrix}1 & 0 \\0 & 1\end{bmatrix} + {\sum\limits_{i}^{\;}{\frac{g_{i}({nT})}{C} \cdot \begin{bmatrix}h_{11}^{i} & 0 \\h_{21}^{i} & 0\end{bmatrix}}}}} \\{= \begin{bmatrix}{1 + {\sum\limits_{i}^{\;}{h_{11}^{i}\frac{g_{i}({nT})}{C}}}} & 0 \\{\sum\limits_{i}^{\;}{h_{21}^{i}\frac{g_{i}({nT})}{C}}} & 1\end{bmatrix}}\end{matrix} & (107)\end{matrix}$

and that its inverse is

$\begin{matrix}{{{\left\lbrack {I + {\sum\limits_{i}^{\;}{\frac{g_{i}({nT})}{C} \cdot H_{i} \cdot B}}} \right\rbrack^{- 1} = {\frac{1}{\Delta}\begin{bmatrix}1 & 0 \\{- {\sum\limits_{i}^{\;}{h_{21}^{i}\frac{g_{i}({nT})}{C}}}} & \Delta\end{bmatrix}}},{where}}{\Delta = {1 + {\sum\limits_{i}^{\;}{h_{11}^{i}\frac{g_{i}({nT})}{C}}}}}} & (108)\end{matrix}$

The computation of the state update may proceed as follows (wheree^(A·T)≡G):

-   -   I. Pre-calculate with the desired quantization and store all        coefficients (G, h_(i), and q).    -   II. For each channel, calculate h₁₁ ^(i)g_(i)(nT)/C and h₂₁        ^(i)g_(i)(nT)/C    -   III. Calculate the sums and scaled-by-E_(i) sums

${\sum\limits_{i}^{\;}\frac{E_{i}h_{11}^{i}{g_{i}({nT})}}{C}},{\sum\limits_{i}^{\;}\frac{h_{11}^{i}{g_{i}({nT})}}{C}},{\sum\limits_{i}^{\;}\frac{E_{i}h_{21}^{i}{g_{i}({nT})}}{C}},{\sum\limits_{i}^{\;}\frac{h_{21}^{i}{g_{i}({nT})}}{C}}$

-   -   IV. Evaluate the right-hand side of Eq. (105) and store in        auxiliary variables, i.e.,

${\begin{bmatrix}g_{11} & g_{12} \\g_{21} & g_{22}\end{bmatrix}\begin{bmatrix}{v({nT})} \\{u({nT})}\end{bmatrix}} + \begin{bmatrix}{\sum\limits_{i}^{\;}\frac{E_{i}h_{11}^{i}{g_{i}({nT})}}{C}} \\{\sum\limits_{i}^{\;}\frac{E_{i}h_{21}^{i}{g_{i}({nT})}}{C}}\end{bmatrix} + \begin{bmatrix}q_{1} \\q_{2}\end{bmatrix}$${auxv} = {{g_{11}{v({nT})}} + {g_{12}{u({nT})}} + {\sum\limits_{i}^{\;}\frac{E_{i}h_{11}^{i}{g_{i}({nT})}}{C}} + q_{1}}$${auxu} = {{g_{21}{v({nT})}} + {g_{22}{u({nT})}} + {\sum\limits_{i}^{\;}\frac{E_{i}h_{21}^{i}{g_{i}({nT})}}{C}} + q_{2}}$

-   -   V. Perform the multiplication by the inverse matrix:

$\begin{matrix}{{v\left( {{nT} + T} \right)} = \frac{auxv}{\Delta}} \\{= \frac{auxv}{1 + {\sum\limits_{i}^{\;}\frac{h_{11}^{i}{g_{i}({nT})}}{C}}}}\end{matrix}$ $\begin{matrix}{{u\left( {{nT} + T} \right)} = {{{- \frac{auxv}{\Delta}}{\sum\limits_{i}^{\;}\frac{h_{21}^{i}{g_{i}({nT})}}{C}}} + {auxu}}} \\{= {{{- {v\left( {{nT} + T} \right)}}{\sum\limits_{i}^{\;}\frac{h_{21}^{i}{g_{i}({nT})}}{C}}} + {auxu}}}\end{matrix}$

Note that for the NMDA case with a double exponential, each exponentialmay most likely be dealt with individually. This leads topre-calculating and storing two vectors, one for the decayingexponential term and one for the rising exponential term.

FIG. 21 is a block diagram 2100 illustrating the realization of Eq.(106). For simplicity, the matrix inversion block 2102 representing thefirst portion of Eq. (106) is not expanded. From the current statevector x(nT), the membrane potential v(nT) from the state vector isquantized in the quantization block 2104 to determine the memory index.The memory index is used to select the predetermined parameters (e.g.,coefficients) for the piecewise linear approximation of the neuron modelcorresponding to the quantization interval of the membrane potential andstored in the memory lookup table(s) 2106. Based on the memory index,the parameters associated with the quantization interval are loaded,including the coefficients for the exponential matrix G 2108 and thedefined vector q 2110 from Eq. (99) above. The parameters also includethe coefficients for input current vector h 2112 and N synaptic currentinput vectors h_(i) 2114. Summation block 2116 sums the input currentI(nT)/C applied to the vector h 2112 and conductance values scaled bythe reverse potentials E_(i)g_(i)(nT)/C applied to the N synapticcurrent input vectors h_(i) 2114. The result is processed according toEq. (106) with the matrix inversion block 2102 to generate the updatedstate vector x(nT+T). A delay block 2118 may be added so that the statesare updated at each step size T in time, rather than continuously.

Other approximations to the integrals in Eq. (101) could be used. Forexample, rather than calculating the exact values and storing them inthe h_(i) vectors, a trapezoidal rule applied to g(t) can be used, thusavoiding storage of the h_(i) vectors and further reducing the memorydemands at the expense of a slightly worse approximation.

SUMMARY

In the subsections above, a general approach to piecewise linearizationof neuron models was described. FIG. 7 summarizes the various steps andlinearization methods. Starting with a nonlinear time varying system at702, piecewise linearization may be performed over the intervalT_(n)≦t≦T_(n+1) at 704 to generate a piecewise linear time-varyingsystem in terms of the matrix A(v_(n),t) at 706. If constantcoefficients (i.e., a piecewise constant function) are used for eachinterval T_(n)≦t≦T_(n+1) at 708, then a piecewise LTI system is formedat 710 in terms of the matrix A(v_(n),T_(n)). If this linearization isperformed independent of the time interval at 712, then a piecewise LTIsystem is created at 714 in terms of the matrix A(v_(n)).

Example State-Space Region-Dependent Linear Neuron Modeling

Certain aspects of the present disclosure described above mainly coverthe partitioning of a two-dimensional state-space limited to one of thestate variables (e.g., the membrane potential v). As described above,the (v,u) state-space is partitioned into vertical blocks using a singlestate, and each block, representing a voltage interval, is associatedwith a distinct set of linear differential equations.

Certain aspects of the present disclosure as described below extend thisconcept. The state space may be partitioned into different regions (notlimited to vertical blocks), where each region is associated with adistinct set of linear differential equations. The regions may benon-overlapping. The description below is primarily focused on thetwo-state (or, equivalently, two-dimensional) model with theunderstanding that the two-state model may be readily extended to threeor more states where the partitions may be volumes and multi-dimensionalregions, respectively.

Review of the Piecewise Linear Formulation as a Function of v

The two-dimensional piecewise linear differential equations may beexpressed as a set of K pairs of differential equations where the n^(th)pair is defined by

$\begin{matrix}{{{C\frac{\;}{t}v} = {{m_{k}v} + d_{k} - u + I}},{v_{k} \leq v < v_{k + 1}}} & (109) \\{{\frac{\;}{t}u} = {{abv} - {au} - {abv}_{r}}} & (110)\end{matrix}$

where kε[1,K] and F(v)=m_(k) v+d_(k). In this approach, the voltagespace is partitioned into K distinct regions, and the k^(th) region isassociated with the k^(th) pair of differential equations. The union ofthe K regions covers the space of (v,u) values. FIG. 8 illustrates thepartitioning given K=4, where the voltage space is divided into fourregions: v<v₁, v₁≦v<v₂, v₂≦v<v₃, and v₃≦v. For K=2, this reduces to theHunzinger Cold neuron model (with the closed and open intervals switchedsuch that v_(k)<v≦v_(k+1)).

Generalized Formulation

The basic idea is to extend the piecewise linear concept beyond regionsdefined solely in terms of v (i.e., allow Eq. (110) to vary with eachinterval, as well as Eq. (109)). This can be accomplished by defining Ktwo-dimensional regions Ω_(k) in the (v,u) space where kε[1,K] and wherethe union of the K regions covers the entire space of (v,u) values.Furthermore, the parameters a and b, in addition to m_(k) and d_(k), areallowed to be distinct for each region such that these parameters becomea function of k. The parameter C may also be made distinct per region asC_(k), but is left as fixed per region henceforth, with theunderstanding that C can be made region-specific. These extensions maynot necessarily have biological equivalents. However, from anengineering and computational perspective, such extensions enable thedesign of richer neuron models and a simpler platform to develop neuronmodels.

The two aforementioned extensions yield:

$\begin{matrix}{{{C\frac{\;}{t}v} = {{m_{k}v} + d_{k} - u + I}},{\left( {v,u} \right) \in \Omega_{k}}} & (111) \\{{{\frac{\;}{t}u} = {{a_{k}b_{k}v} - {a_{k}u} - {a_{k}b_{k}v_{r}}}},{\left( {v,u} \right) \in \Omega_{k}}} & (112)\end{matrix}$

The above expressions may be expressed more compactly, using thestate-space notation described above, as

$\begin{matrix}{{{\overset{.}{x} = {{A_{k}x} + {b\left( \frac{I}{C} \right)} + p_{k}}},{x \in \Omega_{k}}}{where}} & (113) \\{x \equiv \begin{bmatrix}v \\u\end{bmatrix}} & (114) \\{A_{k} \equiv \begin{bmatrix}\frac{m_{k}}{C} & {- \frac{1}{C}} \\{a_{k}b_{k}} & {- a_{k}}\end{bmatrix}} & (115) \\{p_{k} \equiv \begin{bmatrix}{d_{k}/C} \\{{- a_{k}}b_{k}v_{r}}\end{bmatrix}} & (116) \\{b \equiv \begin{bmatrix}1 \\0\end{bmatrix}} & (117)\end{matrix}$

Next, the M-dimensional state-space region-dependent case and thestate-space solution are described.

M-Dimensional State-Space Case

In the M-dimensional case where there are M state variables, thestate-space equations become

$\begin{matrix}{{{\overset{.}{x} = {{A_{k}x} + {b\left( \frac{I}{C} \right)} + p_{k}}},{x \in \Omega_{k}}}{where}} & (118) \\{x \equiv \begin{bmatrix}x_{1} \\\vdots \\x_{M}\end{bmatrix}} & (119) \\{A_{k} \equiv \begin{bmatrix}a_{11}^{(k)} & \ldots & a_{1M}^{(k)} \\\vdots & \ddots & \vdots \\a_{M\; 1}^{(k)} & \ldots & a_{MM}^{(k)}\end{bmatrix}} & (120) \\{p_{k} \equiv \begin{bmatrix}d_{1}^{(k)} \\\vdots \\d_{M}^{(k)}\end{bmatrix}} & (121) \\{b \equiv \begin{bmatrix}1 \\0 \\\vdots \\0\end{bmatrix}} & (122)\end{matrix}$

Region-Dependent State-Space Solutions

The state-space solutions of Eq. (113) in the continuous-time domain anddiscrete-time domain are presented. Their derivations follow naturallyusing the same approach as that described above. The continuous-timesolution may be expressed as

$\begin{matrix}{{{x(t)} = {{^{A_{k}t}{x(0)}} + {\int\limits_{0}^{t}{{^{A_{k}{({t - \tau})}} \cdot b \cdot {{I(\tau)}/C}}{\tau}}} + q_{k}}},{{x(0)} \in {\Omega - k}}} & (123)\end{matrix}$

where x(0) represents the initial value of x(t),

q _(k) ≡A _(k) ⁻¹(e ^(A) ^(k) ^(t) −I)·p _(k)  (124)

and e^(x) represents the matrix exponential of a matrix X.

The corresponding discrete-time solution, with a sampling period of T,may be expressed as:

$\begin{matrix}{{{{x\left( {{nT} + T} \right)} = {{^{A_{k}t} \cdot {x({nT})}} + {h_{k} \cdot {{I({nT})}/C}} + q_{k}}},{{x({nT})} \in \Omega_{k}}}{where}} & (125) \\{h_{k} \equiv {\int\limits_{0}^{T}{^{A_{k}\tau}{\tau}\; b}}} & \left( {126a} \right) \\{= {A_{k}^{- 1}\left( {^{A_{k}t} - 1} \right)}} & \left( {126b} \right)\end{matrix}$

The present state-space vector x(nT) identifies the region Ω_(k) andwhich set of {e^(A) ^(k) ^(t),h_(k),q_(k)} matrix and vectors to use tocompute x(nT+T) according to Eq. (120).

Furthermore, the same solutions for non-synaptic and synaptic currenttypes described above readily apply. This includes the linear-timevarying (LTV) solutions. This would involve replacing thevoltage-dependent regions (vertical partitions) with these more generalstate-space regions.

Illustrative Examples

FIG. 9 illustrates an example with three regions (i.e., Ω_(k) wherekε[1,3]). All three regions in FIG. 9 are rectangular regions.

FIG. 10 illustrates an example with four regions illustrating how theregions need not be rectangular. This particular design may be used tocreate a subthreshold oscillation neuron where Ω₁, Ω₂, Ω₃, and Ω₄represent regions where (v,u), respectively, spirals outward with agrowing radius, spirals in a circular fashion (with no decay or nogrowth), spirals inward with a decaying radius, and tends to spike(equivalent to the ALIF region of the Hunzinger Cold model). This designis described in more detail below.

FIG. 11 illustrates a more academic example with five regions toillustrate that a region with any shape (e.g., a hexagonal, triangular,or star-shaped region) could be created in theory.

Realization of the State-Space Region-Dependent Linear Neuron Modeling

Presently, the regions are defined by

Ω_(k)={(v,u)|v _(k) ≦v<v _(k+1)}

where v_(k) and v_(k+1) represent the left and right boundaries forvoltage. So, at each time step, the index k is determined by taking thepresent state (v(nT),u(nT)) and checking it against each of the Kregions: Ω_(k) kε[1,K]. The index k, represented by k_(F), may becomputed directly via a function of only v as well as the quantizationinterval, the minimum possible voltage, and number of intervals:

$k = {\left\lfloor \frac{{v({nT})} - v_{\min,F}}{\Delta \; v_{F}} \right\rfloor + 1}$

For certain aspects, the regions (or at least a portion thereof) may beslightly overlapping. In this case, hysteresis may be used inidentifying which of the K regions (Ω_(k), kε[1,K]) contains the presentstate (v(nT),u(nT)).

Two examples of generalized regions are provided below. Other regionsmay also be developed.

For Rectangular Regions Dependent on Both v and u

For generalized rectangular regions, the k^(th) region may be defined by

Ω_(k)={(v,u)|v _(k) ≦v<v _(k+1) ,u _(k) ≦u<u _(k+1)}

where v_(k) and v_(k+1) represent the left and right boundaries forvoltage, respectively, and where u_(k) and u_(k+1) represent the lowerand upper boundaries for the recovery variable, respectively. So, ateach time-step, the index k may be determined by identifying which ofthe K regions (Ω_(k), kε[1,K]) contains the present state (v(nT),u(nT)).

For Elliptical Regions

For ease of explanation, assume only elliptical regions each centeredabout the rest state (v_(r),0). Then the k^(th) region may be defined by

Ω_(k)={(v,u)|ρ_(k)≦ρ(v,u)<ρ_(k+1)}

where ρ_(k) is the “radial distance” to the inner elliptical boundaryalong the v-axis, ρ_(k) is the “radial distance” to the outer ellipticalboundary along the v-axis, the “radial distance” may be defined as

ρ(v,u)≡√{square root over ((v−v _(r))²+(u/β)²)}{square root over ((v−v_(r))²+(u/β)²)}

and β is a parameter associated with the ellipse (β=1 leads toconcentric regions). As before, at each time-step, the index k isdetermined by taking the present state (v(nT), u(nT)) and checking itagainst each of the K regions: Ω_(k), kε[1,K].

Advantages of State-Space Region-Dependent Linear Neuron Modeling

Such generalized partitioning of the state-space for neuron modeling canhelp support the synthesis of neurons and do so in a systematic manner.Two examples of this are in implementing: (1) onset detection and eventcounting and (2) subthreshold oscillation behavior, the latter of whichis described in the next section.

An artificial neuron based on the Hunzinger Cold model (which has twopartitions) may be developed for onset detection and event counting.This artificial neuron, initially at rest, is designed to fireimmediately if it receives a spike for onset detection. Thereafter, theneuron is designed to start counting incoming spikes and fire after acertain number of spikes. This continues to repeat itself as spikesarrive. A potential drawback of this design occurs after there are nomore spike arrivals since the neuron can take about 500 ms beforereaching the equilibrium rest state. This may be too long a time beforethe neuron can perform onset detection. Part of the constraint in thisdesign is how shortening the time to the rest state impacts thespike-counting ability by making it behave more as a LIF neuron.

According to certain aspects of the present disclosure, the constraintcan be removed by separating the two behaviors via partitioning,analogous to how the LIF and ALIF behavior of the Izhikevich neuronmodel may be broken by partitioning the two regions with the HunzingerCold model. The time to reset the neuron to perform onset detection maybe reduced by noticing how the u state variable behaves above and belowa threshold of roughly 0.3 in the LIF region. Below that value, it takesa long time for the system to reach the rest state. This can beexpedited by partitioning the state-space in terms of u as well as v asillustrated in the three-region state-space of FIG. 9 to realize animproved onset and event counting neuron. In the LIF region wherev<v_(t), the artificial neuron now behaves as an integrate and fire (IF)neuron for event counting if u>u₀ (i.e., region Ω₂) and a LIF neuronwhich goes to v_(r) otherwise (i.e., region Ω₁). Once the state entersregion Ω₁, the system may be designed to move to the rest state asquickly as desired with shorter time constants. By partitioning the LIFregion as a function of u into the two desired regions (Ω₂ for IF and Ω₁for LIF) with two pairs of differential equations per region, the largedelay may be eliminated.

Example Subthreshold Oscillation Neuron Design Based on State-SpaceRegion-Dependent Linear Neuron Models

Subthreshold oscillations may be used in an artificial nervous system toincrease the likelihood of spiking in the upswing of the membranepotential from the resting potential (i.e., “priming” the artificialneuron to spike) and to decrease the likelihood in the downswing. Theprevious section describes an extension to the piecewise linear neuronmodeling by allowing partitioning of a two-dimensional (ormulti-dimensional) state-space based not only on one of the statevariables (e.g., the voltage), but also a combination and regions (orfunctions) of state variables. This section addresses the shortcomingsof existing subthreshold oscillation neuron designs. For example,existing neuron models can support the design of a subthresholdoscillation neuron exhibiting at most one of three oscillation types(e.g., decaying oscillations as illustrated in the graph 1200 of FIG.12A, sustained oscillations as illustrated in the graph 1210 of FIG.12B, or growing oscillations as illustrated in the graph 1220 of FIG.12C), but these models cannot enable the design of an artificial neuronexhibiting two or more of the three types. Moreover, existing neuronmodels cannot support the systematic design of subthreshold oscillationneurons with a particular oscillation frequency, a rate of decay fordecaying oscillation, a rate of growth for growing oscillations, and/ora magnitude for sustained oscillations.

Accordingly, what is needed are artificial neurons with an improvedneuron model capable of exhibiting all three types of oscillations andsupporting the systematic realization of desired oscillation behaviors.

Comparison Between Subthreshold-Oscillation-Capable Neuron Models

There are four existing neuron models that can produce subthresholdoscillations: the Izhikevich (simple) model, the Adaptive Exponential(AdEx) neuron model (which is the linear-plus-exponential function ofBrette and Gerstner described above), the quartic model, and theHunzinger Cold model. These neuron models may characterized in fourways: (1) which of the three possible subthreshold oscillation behaviors(damped oscillation, sustained oscillation, or growing oscillations)they can produce; (2) how many of the three behaviors can be exhibitedby a single neuron; (3) how easy it is to design a neuron with a desiredoscillation frequency; and (4). how easy it is to fine tune theoscillation design (e.g., defining the magnitude of the sustainedoscillation or defining multiple rates of decays or growths). Thefollowing table summarizes how each neuron model fares for the firstmeasure (where these neuron models can possibly support from one to allof the three oscillation types):

Damped Sustained Growing Neuron Model Oscillations OscillationsOscillations Izhikevich model and Y AdEx model Quartic model Y Y YHunzinger Cold Y Y Y model State-space region- Y Y Y dependent model

The next table summarizes how each neuron model fares for the remainingcharacterization measures:

One Two Three Systematic, Fine type of types of types of simple tuningModel oscillation oscillation oscillation design capability Izhikevich Yand AdEx models Quartic Y model Hunzinger Y Y Cold model State- Y Y Y YY space region- dependent model

Each neuron model (besides the state-space region-dependent approachdescribed below) can reproduce a neuron capable of exhibiting only asingle oscillation type. And only the Hunzinger Cold model, due to itslinear design, is capable of offering a systematic way of designing atarget. None of them (besides the state-space region-dependent approachdescribed below), however, offers the capability of fine tuning.

Subthreshold Oscillation Based on State-Space Region-Dependent LinearNeuron Modeling

Certain aspects of the present disclosure, as indicated in the tablesabove and in contrast with existing neuron models, can: (1) produce anartificial neuron that can exhibit in the subthreshold regime all or anysubset of the three oscillation types and (2) be designed not onlysystematically like the Hunzinger Cold model, but can also support finetuning of the design.

Design

Exhibition of all three oscillation types may be created by using thestate-space region dependent approach described above. For example, afour-region, two-dimensional neuron model may be defined as illustratedin FIG. 10.

Each state-space region Ω_(k) is associated with a set of lineardifferential equations such that

$\begin{matrix}{{{C\frac{\;}{t}v} = {{m_{k}v} + d_{k} - u + I}},{\left( {v,u} \right) \in \Omega_{k}}} & (127) \\{{{\frac{\;}{t}u} = {{a_{k}b_{k}v} - {a_{k}u} - {a_{k}b_{k}v_{r}}}},{\left( {v,u} \right) \in \Omega_{k}}} & (128)\end{matrix}$

or, in state-space vector notation,

$\begin{matrix}{{{\overset{.}{x} = {{A_{k}x} + {b\left( \frac{I}{C} \right)} + p_{k}}},{x \in \Omega_{k}}}{and}} & (129) \\{A_{k} \equiv \begin{bmatrix}\frac{m_{k}}{C} & {- \frac{1}{C}} \\{a_{k}b_{k}} & {- a_{k}}\end{bmatrix}} & (130) \\{p_{k} \equiv \begin{bmatrix}{d_{k}/C} \\{{- a_{k}}b_{k}v_{r}}\end{bmatrix}} & (131)\end{matrix}$

Regions Ω₁ (appearing circular, but is generally elliptical), Ω₂ (alsogenerally elliptical), Ω₃, and Ω₄ may be designed to reproduce growingoscillations, sustained oscillations, decaying oscillations, and ALIFbehavior, respectively.

According to certain aspects, the state-space may be defined by morethan two dimensions. In this case, the following alternate forms may beused instead:

$\begin{matrix}{A_{k} \equiv \begin{bmatrix}a_{11}^{(k)} & \ldots & a_{1M}^{(k)} \\\vdots & \ddots & \vdots \\a_{M\; 1}^{(k)} & \ldots & a_{MM}^{(k)}\end{bmatrix}} & (132) \\{p_{k} \equiv \begin{bmatrix}d_{1}^{(k)} \\\vdots \\d_{M}^{(k)}\end{bmatrix}} & (133)\end{matrix}$

If sustained oscillations are desired, some thickness may be warrantedin case of time-stepped approaches to handle the scenario where v(nT+T)could inadvertently bypass the sustained oscillation region.

Reset Modification

As described above, the reset mechanism may involve resetting u to u+d,which amounts to incrementing u after a spike event or after a controlsignal is activated, for example. For certain aspects, to ensure thatthe u+d value is not so large that it can cause the neuron toimmediately spike or enter sustained oscillation right from the start,it may be desirable to have the option of resetting u to a targetconstant (i.e., to something other than u+d). This way, the designer hascontrol of where in the state-space the neuron resets.

Defining the State-Space Regions

In general, the (v,u) trace in sustained oscillations may take on theform of an ellipse whose major and minor axes may not be aligned withthe v, u Cartesian axes (i.e., the ellipse may have a tilt). Thesolutions may take on the following form:

v(t)=A _(c) ^((v)) cos(wt)+A _(s) ^((v)) sin(wt)  (134a)

u(t)=A _(c) ^((u)) cos(wt)+A _(s) ^((u)) sin(wt)  (134b)

where the four coefficients are functions of the matrix elements in A(where the index k has been dropped) and the initial conditions: (v(0),u(0)).

For ease in exposition and design, the description here focuses oncircular traces whose major and minor axes are equal. These call for thematrix elements to satisfy the following constraints:

$\begin{matrix}\begin{matrix}{a_{11} = a_{22}} \\{= \lambda_{R}} \\{= \frac{1}{\tau}}\end{matrix} & (135) \\\begin{matrix}{a_{12} = {- a_{21}}} \\{= \lambda_{I}} \\{= {2\pi \; f_{c}}}\end{matrix} & (136)\end{matrix}$

where λ_(R) and λ₁ represent the real and imaginary parts, respectively,of the complex eigenvalues: λ_(1,2)=λ_(R)±λ₁, τ represents the (decay orrise) time constant, f_(c) represents the oscillation frequency, and theindex k has been dropped. This results in circular regions and astraightforward form of the matrix:

$\begin{matrix}{A \equiv \begin{bmatrix}\frac{1}{\tau} & {2\pi \; f_{c}} \\{{- 2}\pi \; f_{c}} & \frac{1}{\tau}\end{bmatrix}} & (137)\end{matrix}$

Using this simpler form with circular regions, the radius of the circlecentered about the equilibrium point (v_(r),0) is

ρ(v,u)≡√{square root over ((v−v _(r))² +u ²)}  (138)

representing the distance from (v,u) to the equilibrium point. Forexample, ρ(v_(r)+2 mV,0)=2. So, the k^(th) region may be defined by

Ω_(k)={(v,u)|ρ_(k,in)<ρ(v,u)≦ρ_(k,out) and v≦v _(t)}  (139)

where ρ_(k,in) represents the radius of the inner circular boundary ofΩ_(k) and ρ_(k,out) represents the outer circular boundary of Ω_(k).

For the K=4 example in FIG. 10, the sustained oscillation region Ω₂would be defined by

Ω₂={(v,u)|ρ_(g)<ρ(v,u)≦ρ_(d) and v≦v _(t)}  (140)

where ρ_(g) represents the radius of the circular boundary separating Ω₁and Ω₂ and where ρ_(d) represents that between Ω₂ and Ω₃.

For the growing oscillation region Ω₁:

Ω₁={(v,u)|ρ(v,u)≦ρ_(g) and v≦v _(t)}  (141)

For the decaying oscillation region Ω₃

Ω₃={(v,u)|ρ_(d)<ρ(v,u) and v≦v _(t)}  (142)

The ALIF region Ω₄ in FIG. 10 is defined by

Ω₄={(v,u)|v>v _(t)}  (143)

Additional Regions for Multi-Stage Decay and Growth Regions

The above idea may be extended to support more gradual decaying orgradual growing of the oscillations. For example, to have varying ratesof decaying orbits and varying rates of growing orbits, multiple ringsmay be used as shown in FIG. 13 with six regions. In this example, Ω₃may be considered as the sustained oscillation region. The region Ω₅ mayhave a faster decaying orbit, and Ω₄ may have a slower decaying orbit.This would allow for a graded and more gentle decay to the sustainedoscillation region Ω₃. Similarly, Ω₁ may have a faster growing orbit,and Ω₂ may have a slower growing orbit. This concept may be extended tomore rings for even more gradual changes. Furthermore, multi-stagesustained oscillation regions may be used if different sustainedoscillation parameters (e.g., different frequencies) are desired.

Example Operations for Implementing an Artificial Neuron

FIG. 14 is a flow diagram of example operations 1400 for updating astate of an artificial neuron in an artificial nervous system, inaccordance with certain aspects of the present disclosure. Theoperations 1400 may be performed in hardware (e.g., by one or moreneural processing units (e.g., artificial neurons), such as aneuromorphic processor), in software, or in firmware. The artificialnervous system may be modeled on any of various biological or imaginarynervous systems, such as a visual nervous system, an auditory nervoussystem, the hippocampus, etc.

The operations 1400 may begin, at 1402, by determining that a firststate of the artificial neuron is within a first region. At 1404, asecond state of the artificial neuron is determined based at least inpart on a first set of linear dynamical equations, wherein the first setof linear dynamical equations is based, at least in part, on a first setof parameters (e.g., coefficients) corresponding to the first region. At1406, it is determined that the second state of the artificial neuron iswithin a second region. At 1408, a third state of the artificial neuronis determined based, at least in part, on a second set of lineardynamical equations. The second set of linear dynamical equations may bebased, at least in part, on a second set of parameters corresponding tothe second region.

According to certain aspects, the first and second sets of linearequations include discrete time solutions of linear time-invariant (LTI)state-space equations. For other aspects, the discrete time solutionsmay be based at least in part on closed-form solutions forcontinuous-time state-space equations.

According to certain aspects, the operations 1400 further includefetching at least one of the first or second set of parameters frommemory. As used herein, memory may refer to any of various suitablemeans for storing data, whether permanently or temporarily, locally orremotely to a processing unit, on-chip or off-chip, including randomaccess memory (RAM), cache memory, registers, latches, flip-flops, andthe like. Such fetching may include fetching the at least one of thefirst or second set of parameters from a memory local to the artificialneuron. For certain aspects, the operations 1400 further includecalculating at least a portion of at least one of the first or secondset of parameters. At least the portion of the at least one of the firstor second set of parameters may be calculated using one or more valuesfetched from memory.

According to certain aspects, at least one of the first or second set ofparameters may be obtained by approximating, with a piecewise linearfunction, at least a portion of a nonlinear function in a neuron modelassociated with the artificial neuron. For certain aspects, thenonlinear function includes a membrane potential (v) multiplied with avoltage-dependent conductance (g(v)). The voltage-dependent conductancemay be approximated with a piecewise constant function. For otheraspects, the nonlinear function includes a voltage-dependent function(F(v)), where v is a membrane potential of the artificial neuron. Thepiecewise linear function may have a slope (m) and an intercept (d) foreach of the first and second regions. For certain aspects, the firstregion has a different width in the piecewise linear function than thesecond region. A width of the first or the second region in thepiecewise linear function may be dependent on the nonlinear function.For certain aspects, the piecewise linear approximation is based atleast in part on at least one of a Taylor expansion method, a firstorder linear interpolation method, an optimal linear interpolationmethod, or an average slope method.

According to certain aspects, at least one of the first or second set ofparameters may be designed at least in part to achieve a particularbehavior in the artificial neuron. For example, the first and/or secondset of parameters may be generated (e.g., manually selected by a neuronmodel designer) to effectively create a particular function (e.g., atunable curve), which approximates the desired behavior. In this manner,the artificial neuron may operate based on an entirely new neuron modelor a modification to an existing neuron model.

The first, second, and third states of the artificial neuron may bedefined by a membrane potential (v) and a recovery current (u). Theseare two examples of artificial neuron state variables. According tocertain aspects, the operations 1400 may further include resetting atleast one of the membrane potential or the recovery current of theartificial neuron based at least in part on a determination that a spikeevent has occurred or will occur. The membrane potential may be reset toa resting potential (v(t)→v_(r)). The recovery current may be reset to asum of a present value of the recovery current and an offset(u(t)→u(t)+d).

According to certain aspects, at least one of the first set or secondset of linear equations is based at least in part on a neuron model forthe artificial neuron. For certain aspects, the neuron model is based atleast in part on at least one of an Izhikevich simple model, anexponential-integrate-and-fire (EIF) model, a FitzHugh-Nagumo model, aquartic model, or an intrinsic conductance model (e.g., a neuron modelexpressed as a membrane potential v multiplied with a voltage-dependentconductance g(v)). For certain aspects, the neuron model comprises atleast two dimensions (i.e., at least two state variables). The neuronmodel may be based at least in part on one or more first-order ordinarydifferential equations (ODEs), which may be linear.

According to certain aspects, a step size in time of the neuron model isbased at least in part on a type of the artificial neuron being modeled.A step size in time of the neuron model may be non-uniform. For certainaspects, the step size in time of a particular time step is determinedbased at least in part on at least one of a current or past state of theartificial neuron or on a particular set of parameters (which may beassociated with the current or past state). For certain aspects, a firststep size in time between determining the first and second states isdifferent than a second step size in time between determining the secondand third states. According to certain aspects, determination of thesecond state of the artificial neuron at 1404 may be performed for afirst time step, and determination of the third state of the artificialneuron at 1408 may be performed for a second time step subsequent to thefirst time step.

According to certain aspects, determining at least one of the second orthird state of the artificial neuron is based at least in part oncurrents input to the artificial neuron. For certain aspects, the inputcurrents comprise at least one of a synaptic current or a general,external current. The synaptic current may be based at least in part ona time-dependent conductance for each of one or more channels and areversal potential for each of the one or more channels. For certainaspects, the synaptic current is based at least in part on apost-synaptic membrane potential of the artificial neuron for anN-methyl-D-aspartate (NMDA) channel. The time-dependent conductance maybe modeled by an exponential function, an alpha function, or adifference-of-exponentials function.

For certain aspects, the first and second regions are the same region.

According to certain aspects, the operations 1400 may further includeoutputting an indication of at least one of the first state, the secondstate, or the third state to a display.

According to certain aspects, at least one of the first region or thesecond region is defined by two or more dimensions. The at least one ofthe first region or the second region is an M-dimensional subspace,where M≧2. The M-dimensional subspace may have any shape or be definedby any boundary (e.g., v>10 mV). For certain aspects, the two or moredimensions are two dimensions, and the M-dimensional subspace has atwo-dimensional shape, such as an ellipse, a circle, a polygon, arectangle, or a square. For other aspects, the two or more dimensionsare three dimensions, and the M-dimensional subspace has athree-dimensional shape, such as a sphere, an ellipsoid, a polygonalprism, a rectangular prism, or a cube. For certain aspects, the firstregion is associated with a different set of linear equations than thesecond region. The first and second regions may be partially overlappingor non-overlapping regions. The first and second regions may have atleast one of different shapes or different sizes. For certain aspects,the first region is associated with a first behavior of the artificialneuron, and the second region is associated with a second behavior ofthe artificial neuron, different from the first behavior. For certainaspects, the first and second regions have layered shapes, such asconcentric rings, tubes, rectangular frames, and the like.

According to certain aspects, the two or more dimensions are composed oftwo dimensions defined by a membrane potential (v) and a recoverycurrent (u). For certain aspects, the operations 1400 may furtherinclude resetting at least one of the membrane potential or the recoverycurrent of the artificial neuron based at least in part on at least oneof receipt of a control signal or a determination that a spike event hasoccurred or will occur. In this case, the membrane potential may bereset to a resting potential, and the recovery current may be reset to aconstant.

According to certain aspects, at least a portion of the second regionoverlaps the first region. In this case, determining that the secondstate of the artificial neuron is within the second region at 1406 maybe based at least in part on hysteresis.

According to certain aspects, the artificial neuron is configured foronset detection and event counting. In this case, the first region maybe associated with an integrate-and-fire (IF) behavior for the eventcounting, and the second region may be associated with a leakyintegrate-and-fire (LIF) behavior.

According to certain aspects, the artificial neuron is configured forsubthreshold oscillation behavior. The subthreshold oscillation behaviormay include at least one of a damped oscillation behavior, a sustainedoscillation behavior, or a growing oscillation behavior. For certainaspects, the first region may be associated with the growing oscillationbehavior, the second region may be associated with the sustainedoscillation behavior, and a third region may be associated with thedecaying oscillation behavior. Furthermore, a fourth region isassociated with an ALIF behavior. For certain aspects, the first regionmay have an elliptical shape, whereas the second region may have a ringshape.

Example Common and Flexible Neural Architecture

Certain aspects of the present disclosure generally relate to developinga common and flexible architecture supporting dynamical neuron models.The design goals include low complexity, accurate modeling of the neurondynamics, and the ability to implement any neuron model. With thisflexible architecture, any of various suitable neuron models may beexecuted and substituted as desired. For example, the neuron model mayinclude at least one of an Izhikevich simple model, anexponential-integrate-and-fire (EIF) model, a FitzHugh-Nagumo model, aquartic model, or a Hunzinger Cold model. Such neuron models may beimplemented using the piecewise linear approximations and parameters(e.g., coefficients) associated therewith, as described above.Furthermore, this flexible architecture allows different artificialneurons in an artificial nervous system to operate with different neuronmodels simultaneously.

FIG. 15A illustrates implementation of such a common and flexible neuralarchitecture for a single neural processing unit 1502, where parametersfor neuron models can be selected, loaded, accessed, added, deleted,adjusted, and/or updated, in accordance with certain aspects of thepresent disclosure. The concepts described herein may easily be expandedto a system of neural processing units 1502 (e.g., an artificial nervoussystem), but for ease of description, only a single neural processingunit is shown in FIG. 15A. The neural processing unit 1502 may implementa state machine 1506 that receives one or more inputs and a set ofparameters for a neuron model and outputs one or more state variables asshown. The state machine 1506 may implement any suitable processingalgorithm for an artificial neuron.

The parameters for a given neuron model may be selected from a set ofparameters 1504 for a plurality of neuron models. For certain aspects,for example, an operator may wish to use a certain neuron model in agiven neural processing unit and select this model from a list ofavailable neuron models. For other aspects, a certain neuron model maybe selected by the system of neural processing units based on learningor training operations. The parameters for a given neuron model may beparameters based on the piecewise linear approximations described above.The set of parameters 1504 may be stored in a memory associated withand/or local to the neural processing unit 1502 for certain aspects. Forother aspects, the set of parameters 1504 may be stored in a memoryglobally accessible by multiple neural processing units or a cachememory internal to the neural processing unit 1502. The parameters forneuron model A may be stored in a first memory location 1508 (e.g., anaddressable memory block), while the parameters for neuron model B maybe stored in a second memory location 1510.

In FIG. 15A, the parameters for neuron model A have been loaded into theneural processing unit 1502. As used herein, the term “loading” may bedefined broadly and may include fetching the parameters for a givenneuron model from a memory accessible by the neural processing unit (orby the system of neural processing units); storing the parameters in amemory local to the neural processing unit; or accessing, in a memory,one or more memory areas associated with (e.g., designated for) theneural processing unit. According to certain aspects, loading theparameters for a certain neuron model into a particular neuralprocessing unit may occur in response to a configuration event. Forexample, the configuration event may include powering up the particularneural processing unit, powering up one or more regions of neuralprocessing units (which may or may not include the particular processingunit), or powering up the entire system of neural processing units.

In FIG. 15B, the parameters for neuron model B have been loaded into theneural processing unit 1502. From then on, the state machine 1506 mayoperate based at least in part on these most recently loaded parameters.In this manner, the neural processing unit 1502 may function accordingto a different neuron model, simply by loading different parameters.Furthermore, the parameters for a particular neuron model may be updatedor deleted at any time.

In FIG. 15C, parameters for neuron model C may be added to the set ofparameters 1504. For example, neuron model C may be a recentlydeveloped, purchased, or licensed neuron model that was not availablewhen the neural processing unit 1502 became operational. In otherexample scenarios, the architect or system designer may have notconsidered neuron model C or may have originally thought this model didnot fit the desired application, but now desires to add this model.

In FIG. 15D, the parameters for neuron model C may be stored in a thirdmemory location 1512. The parameters for neuron model C may be loadedinto the neural processing unit 1502, such that the state machine 1506may operate based at least in part on these most recently loadedparameters. The ability to effectively change neuron models for a neuralprocessing unit by simply loading different parameters therein providesa very flexible architecture that can be updated and changed as desired.

FIG. 16 is a flow diagram of example operations 1600 for implementing acombination of a plurality of neuron models in a system of neuralprocessing units (e.g., an artificial nervous system), in accordancewith certain aspects of the present disclosure. The operations 1600 maybe performed in hardware (e.g., by one or more processing units, such asa neuromorphic processor), in software, or in firmware. The artificialnervous system may be modeled on any of various biological or imaginarynervous systems, such as a visual nervous system, an auditory nervoussystem, the hippocampus, etc.

The plurality of neuron models may include any combination of the neuronmodels described above, as well as any suitable neuron models inexistence or as yet undeveloped or undisclosed. For example, theplurality of neuron models may include at least one of an Izhikevichsimple model, an exponential-integrate-and-fire (EIF) model, aFitzHugh-Nagumo model, a quartic model, a Hunzinger Cold model, or anintrinsic conductance model. As used herein, the phrase “combination ofthe neuron models” generally refers to a set of the plurality of neuronmodels, where the set may include a single member. In other words, thecombination of the plurality of neuron models includes one of theplurality and any combination of any subset of the plurality.

The operations 1600 may begin, at 1602, by loading parameters for afirst neuron model selected from the plurality of neuron models into afirst neural processing unit (e.g., an artificial neuron). Theparameters for a given neuron model may be parameters derived from thepiecewise linear approximations of the neuron models as described above.At 1604, a first state of the first neural processing unit may bedetermined based, at least in part, on the loaded parameters for thefirst neuron model. At 1606, a second state of the first neuralprocessing unit may be determined based, at least in part, on theparameters for the first neuron model and on the first state.

According to certain aspects, the plurality of neuron models may beupdated, at 1608. As described above, updating the plurality of neuronmodels may include deleting or adjusting certain parameters for theexisting neuron models or adding parameters for another neuron model(e.g., a new neuron model that was not available at the time theplurality of neuron models were made available to the system or adesired modification to an existing neuron model to effectively create anew neuron model). For certain aspects, updating the plurality of neuronmodels may occur before loading (or re-loading) the parameters for thefirst neuron model and determining the next state of the first neuralprocessing unit.

At 1610, parameters for a second neuron model may optionally be loadedinto the first neural processing unit. These parameters for the secondneuron model may replace the parameters loaded at 1602 for the firstneuron model. The plurality of neuron models need not be updated at 1608before the parameters for the second neuron model are loaded at 1610. At1612, a third state of the first neural processing unit (e.g.,subsequent to the second state) may be determined based, at least inpart, on the parameters for the second neuron model.

According to certain aspects, loading the parameters for the first modelat 1602 (or for the second neuron model at 1610) is in response to aconfiguration event. For certain aspects, the configuration event is apower up for at least one of the system of neural processing units orthe first neural processing unit.

According to certain aspects, determining the third state at 1612 isfurther based, at least in part, on the second state. If a determinationof a subsequent state based on parameters for a subsequent neuron modelis based on a previous state determined under parameters for a previousneuron model (different from the subsequent neuron model), this may bereferred to as a “hot swap.” For certain aspects, the second neuronmodel is configured, at least in part, to generate a subthresholdoscillation under this “hot swap” condition. This oscillation may startfrom the second state or at another time. For certain aspects, theloading of the parameters for the second neuron model at 1610 is based,at least in part, on at least one of a state condition (e.g., the stateentering a particular region or exhibiting a certain behavior, such as arate of change, etc.), a time condition (e.g., based on a timer), or atrigger. The trigger may be external to the first neural processingunit, for example.

According to certain aspects, the operations 1600 may further includeloading parameters for a second neuron model selected from the pluralityof neuron models into a second neural processing unit. Then, a state ofthe second neural processing unit may be determined based, at least inpart, on the parameters for the second neuron model. For certainaspects, the first and second neural processing units are located indifferent areas of the system. For certain aspects, the first and secondneuron models represent different types of neurons. In other words, thedifferent neuron models in the system of neural processing units may beused to emulate neurons that are in different regions of a neuralsystem. For example, neurons in the visual, auditory, or motor controlsystem may be represented by different neuron models.

Loading the parameters for the first, the second, or any other neuronmodel may involve any of various suitable operations. The descriptionbelow includes loading the parameters for the first neuron model forease of description, but such loading applies to the parameters for anyother neuron model, as well. According to certain aspects, loading theparameters for the first neuron model includes loading the parametersfor the first neuron model into a plurality of neural processing unitsincluding the first neural processing unit. For certain aspects, loadingthe parameters for the first neuron model includes fetching at least aportion of the parameters for the first neuron model from a memory. Thismemory may be accessible by the system of neural processing units forcertain aspects. For certain aspects, at least a portion of the memorymay be local to the first neural processing unit. At least portion ofthe memory may be cache memory. For certain aspects, loading theparameters for the first neuron model further includes accessing, in thememory, a memory block associated with the first neuron model. In thiscase, the memory block may be designated by a pointer to a memoryaddress in the memory block.

FIG. 17 illustrates an example block diagram 1700 of components forimplementing the aforementioned methods for operating an artificialneuron using a general-purpose processor 1702 in accordance with certainaspects of the present disclosure. Variables (neural signals), synapticweights, and/or system parameters associated with a computationalnetwork (neural network) may be stored in a memory block 1704, whileinstructions related executed at the general-purpose processor 1702 maybe loaded from a program memory 1706. In an aspect of the presentdisclosure, the instructions loaded into the general-purpose processor1702 may comprise code for determining that a first state of anartificial neuron is within a first region; code for determining asecond state of the artificial neuron based at least in part on a firstset of linear equations, wherein the first set of linear equations isbased at least in part on a first set of parameters corresponding to thefirst region; code for determining that the second state of theartificial neuron is within a second region; and code for determining athird state of the artificial neuron based at least in part on a secondset of linear equations, wherein the second set of linear equations isbased at least in part on a second set of parameters corresponding tothe second region. For certain aspects, at least one of the first orsecond region is defined by two or more dimensions.

In another aspect of the present disclosure, the instructions loadedinto the general-purpose processor 1702 may comprise code for loadingparameters for a first neuron model selected from a plurality of neuronmodels into a first neural processing unit (in a system of neuralprocessing units) and code for determining a first state of the firstneural processing unit based at least in part on the parameters for thefirst neuron model.

FIG. 18 illustrates an example block diagram 1800 of components forimplementing the aforementioned method for operating an artificialneuron where a memory 1802 can be interfaced via an interconnectionnetwork 1804 with individual (distributed) processing units (neuralprocessors) 1806 of a computational network (neural network) inaccordance with certain aspects of the present disclosure. Variables(neural signals), synaptic weights, and/or system parameters associatedwith the computational network (neural network) may be stored in thememory 1802, and may be loaded from the memory 1802 via connection(s) ofthe interconnection network 1804 into each processing unit (neuralprocessor) 1806. In an aspect of the present disclosure, the processingunit 1806 may be configured to determine that a first state of anartificial neuron is within a first region; to determine a second stateof the artificial neuron based at least in part on a first set of linearequations, wherein the first set of linear equations is based at leastin part on a first set of parameters corresponding to the first region;to determine that the second state of the artificial neuron is within asecond region; and to determine a third state of the artificial neuronbased at least in part on a second set of linear equations, wherein thesecond set of linear equations is based at least in part on a second setof parameters corresponding to the second region. For certain aspects,at least one of the first or second region is defined by two or moredimensions.

In another aspect of the present disclosure, the processing unit 1806may be configured to load parameters for a first neuron model selectedfrom a plurality of neuron models into a first neural processing unit(in a system of neural processing units) and to determine a first stateof the first neural processing unit based at least in part on theparameters for the first neuron model.

FIG. 19 illustrates an example block diagram 1900 of components forimplementing the aforementioned method for operating an artificialneuron based on distributed weight memories 1902 and distributedprocessing units (neural processors) 1904 in accordance with certainaspects of the present disclosure. As illustrated in FIG. 19, one memorybank 1902 may be directly interfaced with one processing unit 1904 of acomputational network (neural network), wherein that memory bank 1902may store variables (neural signals), synaptic weights, and/or systemparameters associated with that processing unit (neural processor) 1904.In an aspect of the present disclosure, the processing unit(s) 1904 maybe configured to determine that a first state of an artificial neuron iswithin a first region; determine a second state of the artificial neuronbased at least in part on a first set of linear equations, wherein thefirst set of linear equations is based at least in part on a first setof parameters corresponding to the first region; determine that thesecond state of the artificial neuron is within a second region; anddetermine a third state of the artificial neuron based at least in parton a second set of linear equations, wherein the second set of linearequations is based at least in part on a second set of parameterscorresponding to the second region. For certain aspects, at least one ofthe first or second region is defined by two or more dimensions.

In another aspect of the present disclosure, the processing unit(s) 1904may be configured to load parameters for a first neuron model selectedfrom a plurality of neuron models into a first neural processing unit(in a system of neural processing units) and to determine a first stateof the first neural processing unit based at least in part on theparameters for the first neuron model.

FIG. 20 illustrates an example implementation of a neural network 2000in accordance with certain aspects of the present disclosure. Asillustrated in FIG. 20, the neural network 2000 may comprise a pluralityof local processing units 2002 that may perform various operations ofmethods described above. Each processing unit 2002 may comprise a localstate memory 2004 and a local parameter memory 2006 that storeparameters of the neural network. In addition, the processing unit 2002may comprise a memory 2008 with a local (neuron) model program, a memory2010 with a local learning program, and a local connection memory 2012.Furthermore, as illustrated in FIG. 20, each local processing unit 2002may be interfaced with a unit 2014 for configuration processing that mayprovide configuration for local memories of the local processing unit,and with routing connection processing elements 2016 that providerouting between the local processing units 2002.

According to certain aspects of the present disclosure, each localprocessing unit 2002 may be configured to determine parameters of theneural network based upon desired one or more functional features of theneural network, and develop the one or more functional features towardsthe desired functional features as the determined parameters are furtheradapted, tuned and updated.

The various operations of methods described above may be performed byany suitable means capable of performing the corresponding functions.The means may include various hardware and/or software component(s)and/or module(s), including, but not limited to a circuit, anapplication specific integrated circuit (ASIC), or processor. Forexample, the various operations may be performed by one or more of thevarious processors shown in FIGS. 17-20. Generally, where there areoperations illustrated in figures, those operations may havecorresponding counterpart means-plus-function components with similarnumbering. For example, operations 1400 illustrated in FIG. 14correspond to means 1400A illustrated in FIG. 14A.

For example, means for displaying may comprise a display (e.g., amonitor, flat screen, touch screen, and the like), a printer, or anyother suitable means for outputting data for visual depiction (e.g., atable, chart, or graph). Means for processing, means for generating,means for loading, means for resetting, means for fetching, means forupdating, means for calculating, means for computing, means foroutputting, or means for determining may comprise a processing system,which may include one or more processors or processing units. Means forstoring may comprise a memory or any other suitable storage device(e.g., RAM), which may be accessed by the processing system.

As used herein, the term “determining” encompasses a wide variety ofactions. For example, “determining” may include calculating, computing,processing, deriving, investigating, looking up (e.g., looking up in atable, a database or another data structure), ascertaining, and thelike. Also, “determining” may include receiving (e.g., receivinginformation), accessing (e.g., accessing data in a memory), and thelike. Also, “determining” may include resolving, selecting, choosing,establishing, and the like.

As used herein, a phrase referring to “at least one of” a list of itemsrefers to any combination of those items, including single members. Asan example, “at least one of a, b, or c” is intended to cover a, b, c,a-b, a-c, b-c, and a-b-c.

The various illustrative logical blocks, modules, and circuits describedin connection with the present disclosure may be implemented orperformed with a general purpose processor, a digital signal processor(DSP), an application specific integrated circuit (ASIC), a fieldprogrammable gate array signal (FPGA) or other programmable logic device(PLD), discrete gate or transistor logic, discrete hardware componentsor any combination thereof designed to perform the functions describedherein. A general-purpose processor may be a microprocessor, but in thealternative, the processor may be any commercially available processor,controller, microcontroller, or state machine. A processor may also beimplemented as a combination of computing devices, e.g., a combinationof a DSP and a microprocessor, a plurality of microprocessors, one ormore microprocessors in conjunction with a DSP core, or any other suchconfiguration.

The steps of a method or algorithm described in connection with thepresent disclosure may be embodied directly in hardware, in a softwaremodule executed by a processor, or in a combination of the two. Asoftware module may reside in any form of storage medium that is knownin the art. Some examples of storage media that may be used includerandom access memory (RAM), read only memory (ROM), flash memory, EPROMmemory, EEPROM memory, registers, a hard disk, a removable disk, aCD-ROM and so forth. A software module may comprise a singleinstruction, or many instructions, and may be distributed over severaldifferent code segments, among different programs, and across multiplestorage media. A storage medium may be coupled to a processor such thatthe processor can read information from, and write information to, thestorage medium. In the alternative, the storage medium may be integralto the processor.

The methods disclosed herein comprise one or more steps or actions forachieving the described method. The method steps and/or actions may beinterchanged with one another without departing from the scope of theclaims. In other words, unless a specific order of steps or actions isspecified, the order and/or use of specific steps and/or actions may bemodified without departing from the scope of the claims.

The functions described may be implemented in hardware, software,firmware, or any combination thereof. If implemented in hardware, anexample hardware configuration may comprise a processing system in adevice. The processing system may be implemented with a busarchitecture. The bus may include any number of interconnecting busesand bridges depending on the specific application of the processingsystem and the overall design constraints. The bus may link togethervarious circuits including a processor, machine-readable media, and abus interface. The bus interface may be used to connect a networkadapter, among other things, to the processing system via the bus. Thenetwork adapter may be used to implement signal processing functions.For certain aspects, a user interface (e.g., keypad, display, mouse,joystick, etc.) may also be connected to the bus. The bus may also linkvarious other circuits such as timing sources, peripherals, voltageregulators, power management circuits, and the like, which are wellknown in the art, and therefore, will not be described any further.

The processor may be responsible for managing the bus and generalprocessing, including the execution of software stored on themachine-readable media. The processor may be implemented with one ormore general-purpose and/or special-purpose processors. Examples includemicroprocessors, microcontrollers, DSP processors, and other circuitrythat can execute software. Software shall be construed broadly to meaninstructions, data, or any combination thereof, whether referred to assoftware, firmware, middleware, microcode, hardware descriptionlanguage, or otherwise. Machine-readable media may include, by way ofexample, RAM (Random Access Memory), flash memory, ROM (Read OnlyMemory), PROM (Programmable Read-Only Memory), EPROM (ErasableProgrammable Read-Only Memory), EEPROM (Electrically ErasableProgrammable Read-Only Memory), registers, magnetic disks, opticaldisks, hard drives, or any other suitable storage medium, or anycombination thereof. The machine-readable media may be embodied in acomputer-program product. The computer-program product may comprisepackaging materials.

In a hardware implementation, the machine-readable media may be part ofthe processing system separate from the processor. However, as thoseskilled in the art will readily appreciate, the machine-readable media,or any portion thereof, may be external to the processing system. By wayof example, the machine-readable media may include a transmission line,a carrier wave modulated by data, and/or a computer product separatefrom the device, all which may be accessed by the processor through thebus interface. Alternatively, or in addition, the machine-readablemedia, or any portion thereof, may be integrated into the processor,such as the case may be with cache and/or general register files.

The processing system may be configured as a general-purpose processingsystem with one or more microprocessors providing the processorfunctionality and external memory providing at least a portion of themachine-readable media, all linked together with other supportingcircuitry through an external bus architecture. Alternatively, theprocessing system may be implemented with an ASIC (Application SpecificIntegrated Circuit) with the processor, the bus interface, the userinterface, supporting circuitry, and at least a portion of themachine-readable media integrated into a single chip, or with one ormore FPGAs (Field Programmable Gate Arrays), PLDs (Programmable LogicDevices), controllers, state machines, gated logic, discrete hardwarecomponents, or any other suitable circuitry, or any combination ofcircuits that can perform the various functionality described throughoutthis disclosure. Those skilled in the art will recognize how best toimplement the described functionality for the processing systemdepending on the particular application and the overall designconstraints imposed on the overall system.

The machine-readable media may comprise a number of software modules.The software modules include instructions that, when executed by theprocessor, cause the processing system to perform various functions. Thesoftware modules may include a transmission module and a receivingmodule. Each software module may reside in a single storage device or bedistributed across multiple storage devices. By way of example, asoftware module may be loaded into RAM from a hard drive when atriggering event occurs. During execution of the software module, theprocessor may load some of the instructions into cache to increaseaccess speed. One or more cache lines may then be loaded into a generalregister file for execution by the processor. When referring to thefunctionality of a software module below, it will be understood thatsuch functionality is implemented by the processor when executinginstructions from that software module.

If implemented in software, the functions may be stored or transmittedover as one or more instructions or code on a computer-readable medium.Computer-readable media include both computer storage media andcommunication media including any medium that facilitates transfer of acomputer program from one place to another. A storage medium may be anyavailable medium that can be accessed by a computer. By way of example,and not limitation, such computer-readable media can comprise RAM, ROM,EEPROM, CD-ROM or other optical disk storage, magnetic disk storage orother magnetic storage devices, or any other medium that can be used tocarry or store desired program code in the form of instructions or datastructures and that can be accessed by a computer. Also, any connectionis properly termed a computer-readable medium. For example, if thesoftware is transmitted from a website, server, or other remote sourceusing a coaxial cable, fiber optic cable, twisted pair, digitalsubscriber line (DSL), or wireless technologies such as infrared (IR),radio, and microwave, then the coaxial cable, fiber optic cable, twistedpair, DSL, or wireless technologies such as infrared, radio, andmicrowave are included in the definition of medium. Disk and disc, asused herein, include compact disc (CD), laser disc, optical disc,digital versatile disc (DVD), floppy disk, and Blu-ray® disc where disksusually reproduce data magnetically, while discs reproduce dataoptically with lasers. Thus, in some aspects computer-readable media maycomprise non-transitory computer-readable media (e.g., tangible media).In addition, for other aspects computer-readable media may comprisetransitory computer-readable media (e.g., a signal). Combinations of theabove should also be included within the scope of computer-readablemedia.

Thus, certain aspects may comprise a computer program product forperforming the operations presented herein. For example, such a computerprogram product may comprise a computer readable medium havinginstructions stored (and/or encoded) thereon, the instructions beingexecutable by one or more processors to perform the operations describedherein. For certain aspects, the computer program product may includepackaging material.

Further, it should be appreciated that modules and/or other appropriatemeans for performing the methods and techniques described herein can bedownloaded and/or otherwise obtained by a device as applicable. Forexample, such a device can be coupled to a server to facilitate thetransfer of means for performing the methods described herein.Alternatively, various methods described herein can be provided viastorage means (e.g., RAM, ROM, a physical storage medium such as acompact disc (CD) or floppy disk, etc.), such that a device can obtainthe various methods upon coupling or providing the storage means to thedevice. Moreover, any other suitable technique for providing the methodsand techniques described herein to a device can be utilized.

It is to be understood that the claims are not limited to the preciseconfiguration and components illustrated above. Various modifications,changes and variations may be made in the arrangement, operation anddetails of the methods and apparatus described above without departingfrom the scope of the claims.

What is claimed is:
 1. A method for implementing a combination of aplurality of neuron models in a system of neural processing units,comprising: loading parameters for a first neuron model selected fromthe plurality of neuron models into a first neural processing unit;determining a first state of the first neural processing unit based atleast in part on the parameters for the first neuron model; anddetermining a second state of the first neural processing unit based atleast in part on the parameters for the first neuron model and on thefirst state.
 2. The method of claim 1, further comprising updating theplurality of neuron models.
 3. The method of claim 1, furthercomprising: loading parameters for a second neuron model into the firstneural processing unit; and determining a third state of the firstneural processing unit based at least in part on the parameters for thesecond neuron model.
 4. The method of claim 3, wherein loading theparameters for the second neuron model is in response to a configurationevent.
 5. The method of claim 4, wherein the configuration eventcomprises a power up for at least one of the system of neural processingunits or the first neural processing unit.
 6. The method of claim 3,wherein determining the third state is further based at least in part onthe second state.
 7. The method of claim 6, wherein the second neuronmodel is configured at least in part to generate an oscillation.
 8. Themethod of claim 6, wherein loading the parameters for the second neuronmodel is based at least in part on at least one of a state condition, atime condition, or a trigger, wherein the trigger is external to thefirst neural processing unit.
 9. The method of claim 1, furthercomprising: loading parameters for a second neuron model selected fromthe plurality of neuron models into a second neural processing unit; anddetermining a state of the second neural processing unit based at leastin part on the parameters for the second neuron model.
 10. The method ofclaim 9, wherein the first and second neuron models represent differenttypes of neurons.
 11. The method of claim 1, wherein loading theparameters for the first neuron model comprises fetching at least aportion of the parameters for the first neuron model from a memory. 12.The method of claim 11, wherein at least a portion of the memory islocal to the first neural processing unit.
 13. The method of claim 11,wherein at least a portion of the memory is cache memory.
 14. The methodof claim 11, wherein loading the parameters for the first neuron modelfurther comprises accessing, in the memory, a memory block associatedwith the first neuron model.
 15. The method of claim 14, wherein thememory block is designated by a pointer to a memory address in thememory block.
 16. The method of claim 1, wherein the plurality of neuronmodels comprises at least one of an Izhikevich simple model, anexponential-integrate-and-fire (EIF) model, a FitzHugh-Nagumo model, aquartic model, or an intrinsic conductance model.
 17. An apparatus forimplementing a combination of a plurality of neuron models in a systemof neural processing units, comprising: a processing system configuredto: load parameters for a first neuron model selected from the pluralityof neuron models into a first neural processing unit; determine a firststate of the first neural processing unit based at least in part on theparameters for the first neuron model; and determine a second state ofthe first neural processing unit based at least in part on theparameters for the first neuron model and on the first state; and amemory coupled to the processing system.
 18. The apparatus of claim 17,wherein the processing system is further configured to update theplurality of neuron models.
 19. The apparatus of claim 17, wherein theprocessing system is further configured to: load parameters for a secondneuron model into the first neural processing unit; and determine athird state of the first neural processing unit based at least in parton the parameters for the second neuron model.
 20. The apparatus ofclaim 19, wherein the processing system is configured to load theparameters for the second neuron model in response to a configurationevent.
 21. The apparatus of claim 20, wherein the configuration eventcomprises a power up for at least one of the system of neural processingunits or the first neural processing unit.
 22. The apparatus of claim19, wherein the processing system is configured to determine the thirdstate further based at least in part on the second state.
 23. Theapparatus of claim 22, wherein the second neuron model is configured atleast in part to generate an oscillation.
 24. The apparatus of claim 22,wherein the processing system is configured to load the parameters forthe second neuron model based at least in part on at least one of astate condition, a time condition, or a trigger, wherein the trigger isexternal to the first neural processing unit.
 25. The apparatus of claim17, wherein the processing system is further configured to: loadparameters for a second neuron model selected from the plurality ofneuron models into a second neural processing unit; and determine astate of the second neural processing unit based at least in part on theparameters for the second neuron model.
 26. The apparatus of claim 25,wherein the first and second neuron models represent different types ofneurons.
 27. The apparatus of claim 17, wherein the processing system isconfigured to load the parameters for the first neuron model by fetchingat least a portion of the parameters for the first neuron model from thememory.
 28. The apparatus of claim 27, wherein at least a portion of thememory is local to the first neural processing unit.
 29. The apparatusof claim 27, wherein at least a portion of the memory is cache memory.30. The apparatus of claim 27, wherein the processing system isconfigured to load the parameters for the first neuron model byaccessing, in the memory, a memory block associated with the firstneuron model.
 31. The apparatus of claim 30, wherein the memory block isdesignated by a pointer to a memory address in the memory block.
 32. Theapparatus of claim 17, wherein the plurality of neuron models comprisesat least one of an Izhikevich simple model, anexponential-integrate-and-fire (EIF) model, a FitzHugh-Nagumo model, aquartic model, or an intrinsic conductance model.
 33. An apparatus forimplementing a combination of a plurality of neuron models in a systemof neural processing units, comprising: means for loading parameters fora first neuron model selected from the plurality of neuron models into afirst neural processing unit; means for determining a first state of thefirst neural processing unit based at least in part on the parametersfor the first neuron model; and means for determining a second state ofthe first neural processing unit based at least in part on theparameters for the first neuron model and on the first state.
 34. Theapparatus of claim 33, further comprising: means for loading parametersfor a second neuron model into the first neural processing unit; andmeans for determining a third state of the first neural processing unitbased at least in part on the parameters for the second neuron model.35. A computer program product for implementing a combination of aplurality of neuron models in a system of neural processing units,comprising a computer-readable storage device having instructionsexecutable to: load parameters for a first neuron model selected fromthe plurality of neuron models into a first neural processing unit;determine a first state of the first neural processing unit based atleast in part on the parameters for the first neuron model; anddetermine a second state of the first neural processing unit based atleast in part on the parameters for the first neuron model and on thefirst state.